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ABSTRACT 


This thesis investigates linear state space modeling of single main rotor helicopters culminating 
in a computer program that can be used for 1) stability and control analysis for any single main rotor 
helicopter or 2) preliminary design of a helicopter. The trim solution for a flight condition is found, 
the aircraft is perturbed about the nominal point, and the stability and control derivatives arc 
determined. State space models and analysis tools are provided by the program. A notional attack 
helicopter designed for the 1993 American Helicopter Society Design Competition and a notional 
utility helicopter are used as examples. 
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I. INTRODUCTION 


A. NATURE OF THE PROBLEM 

M Life is too short to worry about helicopter dynamics.” 

Louis V. Schmidt 
Professor Emeritus 
Nava] Postgraduate School 

Helicopter dynamics, because of the complexity added by the high-mass, rotational components 
of the aircraft, namely the main and tail rotors, is often shunned by the aerodynamicist. The fact that 
main and tail rotors are the dominant contributors to flight dynamics makes the situation even worse. 
The understanding of flight mechanics, stability, and control is difficult enough without this added 
twist. As a result, courses in flight mechanics focus on solutions for fixed wing aircraft. 

As will be shown, the flight mechanics of a helicopter are essentially no different than those 
of a fixed wing aircraft. The analysis of a helicopter in flight depends on solving the equations of 
motion, just like for a fixed wing aircraft. The key difference is the realization of a trimmed solution 
for the helicopter, which is an iterative process vice a straight forward, closed form solution as for 
the fixed wing aircraft. Once the trim solution is found, the evaluation of stability and control 
derivatives for the helicopter is done the same as for a fixed wing aircraft except that the derivatives 
for the main and tail rotors must also be evaluated. 

An attempt was made in the preparation of this paper to be as consistent with fixed wing 
stability notation as possible. The large number of helicopter parameters which are unique or differ 
from fixed wing notation are listed in Appendix A. Additionally, the symbol R is used for rolling 
moment instead of the more common L to be consistent with Ray Prouty’s notation, common to 
helicopter literature. 
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B. BACKGROUND 


This project is the result of two needs at the Naval Postgraduate School. 

1. MODELING OF ROTORCRAFT 

Modeling of rotorcraft needs to be performed because the school is pursuing the 
acquisition of aircraft for flight testing of an auto-land helicopter system. The key deficiency in the 
use of simple, off-the-shelf rotorcraft is the lack of design data and theoretical analysis available for 
control system design. The modeling presented in this thesis is linear. Non-linear modelling can be 
performed on Flightlab which is available at the school. 

2. SOFTWARE AVAILABILITY 

Software available for preliminary helicopter design, such as HESCOMP or CAMRAD, 
is too complicated for students to use during a one quarter Helicopter Design course. The time 
required for a student to learn to use available software is prohibitive. The need for user friendly 
software, based on easily accessible programs such as MATLAB, culminated in two software 
packages being prepared for preliminary helicopter design. The first package JANRAD Performance 
[Ref. 1] solves the helicopter trim solution and is used to perform aerodynamic analyses. The second 
package comes from this study and is used to perform stability and control analysis, but uses the trim 
solution from JANRAD. Both packages require a working knowledge of MATLAB. MATLAB was 
used as the base program for these software packages because it is readily available and familiar to 
students because of its wide use at Naval Postgraduate School. 

3. ASSUMPTIONS 

The following assumptions were made in the foregoing analysis. 

• A person using the program has a working knowledge of MATLAB. 

• The aircraft trim condition is in level flight with climb angle (y c ) equal to zero and no sideslip. 
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• The main and tail rotors provide the only significant contribution to stability and control at a 
hover. Contributions of all other airframe components can be neglected at a hover. 

• Effects of the separate components of the helicopter contribute linearly to stability and control. 

• A single main rotor helicopter is used with a counter-clockwise (as viewed from above) 
rotating main rotor system. 

• The main rotor system is articulated, bearingless or rigid in plane. 

• Rotor blades are rigid. 

• The tail rotor rotates counter-clockwise as viewed from the left side of the helicopter. 

• The tail rotor is a teetering rotor with a delta-3 (pitch-flap coupling) hinge. 

• The tail rotor mast is perpendicular to the direction of flight. 

• Tail rotor coning can be approximated by using only tail rotor thrust and centrifugal forces. 

• The NOTAR” boom is circular. 

• Inflow is uniform. 

• The fuselage of the helicopter is a rigid body and has aerodynamic characteristics similar to 
the Prouty example helicopter [Ref. 2]. 
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n. SOFTWARE OVERVIEW 

A. SOFTWARE ARCHITECTURE 

The Joint Anny/Navy Rotorcraft Analysis and Design program, or JANRAD, consists of two 
major subroutines with others under development. The first routine, JANRAD Performance, 
calculates the trim solution and various performance parameters of a helicopter. It is described in 
detail in Ref. 1. 

The second major component is JANRAD Stability, written by the author, which 1) calculates 
all the stability derivatives of each component of the helicopter, 2) determines the linear state space 
mode] at any trimmable point in the forward flight regime, 3) provides the eigenvalues of the plant 
matrix and plots them in the Argand plane, and 4) determines the open loop transfer functions from 
control inputs to state outputs of the aircraft. 

The basic architecture of the program is in Figure 2.1. The opening menu has the user enter 
data for performance calculations, then branches to either the performance or stability and control 
routines as desired. When the user selects the stability routines, additional input data is requested by 
the program for analysis. Data required for each major aircraft component is delineated under the 
chapters describing the stability and control solution for the component, and is listed in Appendix F. 

Upon data entry, the program performs either a hovering or forward flight stability analysis 
of the helicopter. For calculational purposes, hovering flight is any velocity less than 20 ft/sec or 12 
knots. Each analysis requires a number of perturbations about the nominal flight condition, which 
is done by performing a perturbation, then by calling the TRIM subroutine to determine the new trim 
condition. This procedure is described in more detail in following sections. 


4 







notar? — NOTAR Rotor 1 



Figure 2.1. Program Architecture. 
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After the stability and control analysis is performed, the user can select from five different 
screen outputs depending on his/her needs. Information not available from the output screens can be 
called from the workspace as described in section III.A. 

The linear models provided by the software are saved in a MATLAB data file and can be used 
to perform simulations of aircraft response due to particular control inputs. This is described in the 
Software Use, and Examples chapters. 

B. MODELLING CONVENTIONS 

The coordinate system used in the modelling of the helicopter is the standard set of stability 
and control axes (Fig. 2.2). The x axis passes through the longitudinal axis of the helicopter with 
the positive direction to the front of the aircraft. The y axis passes through the lateral axis of the 
helicopter with the positive direction to the right of the aircraft. The z axis passes through the 
vertical axis with the positive direction to the bottom of the aircraft, or toward the earth if the aircraft 
is upright. This set of coordinates is standard in flight mechanics, though it may be a bit confusing 
initially. 

The coordinate system describing the locations of various aircraft components is the same 
coordinate system used to compute weight and balance for an aircraft. The fuselage station (x axis) 
determines the longitudinal position of the component, the buttline (y axis) determines the lateral 
position, and the waterline (z axis) determines the vertical component. The axes are oriented with 
the x axis parallel to the longitudinal axis of the aircraft with positive direction to the rear, the y axis 
is parallel to the lateral axis with the positive direction to the right side of the aircraft. The z 
direction is orthogonal to the x and y axis with the positive direction upward (Fig. 2.2). 

Modelling of the helicopter was done in standard state-space format, where the dynamics of 
the aircraft are described by a system of differential equations, written as the following system of 
matrix equations. 
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X = [A] X * [B] u 
y = [Q X ♦ [D] ii 


X is the state vector which contains the current state of the system. States are selected which 
are the outputs of integrators and that describe the motion of the aircraft. For the complete coupled 
system, the following states were selected per Ref. 3. 

u - velocity of the aircraft in the positive x direction 
w - velocity of the aircraft in the positive z direction 
q - pitch rate of the aircraft 
0 - pitch angle of the aircraft 

v - velocity of the aircraft in the positive y direction 
p - roll rate of the aircraft 
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♦ - roll angle of the aircraft 
r - yaw rate of the aircraft 
¥ - yaw angle of the aircraft 

The first four states are used to describe the longitudinal characteristics of the helicopter while the 
last five are used to describe the lateral-directional characteristics. A fully coupled system is also 
described using these states. 

The vector u contains the inputs to the system. In the case of the helicopter, the inputs are 

5B1 - change in longitudinal cyclic or elevator 

50 „„ - change in collective 

SA1 - change in lateral cyclic or aileron 

dp - change in anti-torque pedals or rudder 

[A] is the plant matrix and describes the dynamics of the aircraft. The plant matrix is 
determined by using Newton's laws to sum forces and moments about the helicopter at an equilibrium 
condition. This yields a system of differential equations. The system of equations is written in matrix 
form where the matrix relating the states X with the derivatives of the states is the plant matrix [A]. 
The plant matrix used is described in NASA Technical Memorandum 84281 [Ref. 3]. 

The [B] matrix shows the contributions of the inputs u to the system. [B] multiplies the input 
vector u, and the result is added to the system. Each element of [B] is the partial derivative of the 
corresponding state with respect to the input. For example B(1,2) would be the partial derivative of 
u with respect to collective, du/30 oM . 

The vector y is the output of the system. The desired output is described by using the [C] and 
[D] matrices. The matrix [C] is used to determine the output of the system from the states X. If the 
values of the state vector are the desired output, as in the case of this paper, then the matrix is an 
identity matrix the size of the plant matrix; however, the [C] matrix does not have to be the identity 
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matrix and can be set up to provide a wide variety of output. For example, if the desired output is 
to be lateral acceleration, or v, then the fifth row of the plant matrix, which corresponds to the 
equation 

dv du dv dw dv dd dv dq dvdv dv dp dv d& dvdr dv dt 

du dt dw dt 86 dt dq dt dv dt dp di 3$ dt dr dt dif dt 

would be used in the [CJ matrix to provide the output. 

The final matrix, [D], is the direct transition matrix for the input. It can be used to provide 
a direct feed from the input vector u to the output vector y. 

C. THE TRIM SOLUTION 

Before a stability and control analysis can be performed on a helicopter, the trim solution for 
the main and tail rotors must be obtained. The trim solution solves for the collective pitch (0,*,), the 
thrust moment of the blade (Tm), the first lateral harmonic of blade feathering (A,), the first 
longitudinal harmonic of blade feathering (B,), the rotor coning angle (ao M ) and the tip path plane 
angle (a^ M ) of the main rotor. When these values are solved for, the first coefficient of longitudinal 
blade flapping (a UM ), and the first coefficient of lateral blade flapping (b lfM ), are then computed. 
Additionally, the trim solution for the tail rotor is needed to provide tail rotor collective (0 oT ), the first 
coefficient of longitudinal blade flapping (a,, T ), the first coefficient of lateral blade flapping (b UT ), and 
the tail rotor coning angle (ao T ) for stability analysis. 

The trim solution for the main rotor is solved by the TRIM routine in JANRAD. The values 
a ltM , t>iiM arc then computed in the Stability and Control Routine. The tail rotor solution is 
approximated using closed form solutions in the Stability and Control Routine. 

D. GENERAL PROCEDURE TO EVALUATE STABILITY DERIVATIVES 

Stability derivatives are determined by using closed form solutions whenever possible, and by 
solving multiple trim solutions about a nominal point to solve for unknowns. 




By using the TRIM routine in die JANRAD program, the trim solution for the desired flight 
condition is determined. The helicopter is then perturbed about this nominal condition and 
retrimmed. From these additional solutions, the basic stability derivatives for the main and tail rotors 
are solved for. 

By solving for the basic stability derivatives first, the helicopter can be described by using the 
same states as in the a standard linear model of a fixed wing aircraft [u v w p q r ♦ fl ♦]'. Without 
solving for the intermediate derivatives first, the helicopter would have to be described using 
additional states for the main rotor, [a liM b„ M C T /cr M C H /o M C Q /a M X M ']', and the tail rotor, [a UT 
b UT C|/o T C H /o T Cq/o t hj X t ]'. Aside from the fact that a model described in this fashion is 
inconsistent with the standard flight mechanic description, the size of the system model becomes large 
and intimidating. This is part of the reason that helicopter flight mechanics are generally avoided. 

After the basic stability derivatives are determined, intermediate derivatives are evaluated. 
Closed form solutions are used whenever possible to solve intermediate derivatives. After an 
intermediate derivative is determined, the basic derivative and the intermediate derivative are 
combined using the chain rule to form the complete stability derivative. An example solution of a 
main rotor derivative follows. 


Wanted : | —| 

30« 


'u 


Find from rotor plots or numerically: 


Compute: 

«) 


dCJo\ 

dCj/a 

da, 

>. 

36 0 ’ 

ae 0 * 

38 0 


da. 




■§] = ~pA b (QR)* 


dCg/a dCg/o da l dC^o 

‘ - 


ae 0 3a. 36 0 


38, 
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m. STABILITY AND CONTROL DERIVATIVES 


A. MAIN ROTOR CONTRIBUTIONS 

The main rotor is the key contributor to helicopter flight dynamics. The strong coupling 
between the states of the rotor produce complicated responses in many of the aircraft states at the 
same time. This can be visualized through the following example. 

If collective pitch is increased in forward flight, the increased drag in the rotor system 
increases the amount of torque required from the engine; this results in a yaw of the helicopter. 
Simultaneously, since the tip path plane is tilted forward, the forward thrust increases which increases 
the advanced ratio. Tl increase in advance ratio results in differing amounts of lift on the advancing 
and retreating sides of the rotor disk which results in a pitching moment because of the phase lag of 
the blades. At the same time, the pitching moment changes the angle of attack on the fore and aft 
blades, which changes lateral flapping, which changes the rolling moment. Because of the cross 
coupling between lateral and longitudinal flapping, there is a change in longitudinal flapping. And 
so on. 

The basic stability derivatives of the main rotor are described in terms of the partial derivatives 
of the standard helicopter coefficients, a„ M , b„ M , Cj/a M , C H /o M , Cg/a M , fi M , X', A IM , B )M , and the 
main rotor blade stiffness 


(dAT | , ( dR \ 

I da. ' rib. 

\ *«/j# V *»/|f 


--A h pK&R?a 
4 R * 

Y 
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These basic main rotor derivatives are then combined as described in die prior section to form the 
main rotor stability derivatives. A detailed list of the basic main rotor derivatives and die stability 
derivatives is available in Appendix C. 

To solve for the derivatives, JAN RAD uses input data from the performance menu, and output 
computations from the TRIM section of the program. The input required for the trim solution is 
found in the JANRAD Performance thesis [Ref. 1]. Output from the performance routine and input 
from the Stability and Control Routine New File menu (Fig. 3.1) are used to solve for a liM , b l(M 
Cf/o M , C H /o M , Cq/o m , Hh, X', A,, B,. Input data can be changed in the Stability and Control 
Additional Parameters Menu screen 1 (Fig. 3.2). 

Blade flapping moment of inertia (slug ft A 2): 

Hub height above reference datum/waterline (ft): 

Hub fuselage station (ft): 

Hub position right of buttline (ft): 

Mast incidence (negative forward - deg): 

Figure 3.1. Main Rotor New File Menu. 

The perturbation solution in conjunction with the equations in Appendix C solve for the partial 
derivatives of a UM , b„ M C^la^, C H /o M , Cq/o m , X', A,, B, with respect to each other, thus 
completing the evaluation of the basic main rotor derivatives. The CRUISE.M and HOVER.M 
routines then combine the basic and intermediate derivatives to give a produce a complete set of main 
rotor derivatives. 

Main rotor output is displayed in output screen 2 (Fig. 3.3). Derivatives are stored in the 
workspace. Because of the large number of derivatives, it is impractical to display all of them on the 
screen. Instead, key derivatives and control power are displayed in the Key Control Parameters 
output screen (Fig. 3.4) AFTER all contributions from other aircraft components are determined. 
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*** STABILITY AND CONTROL MEND *** 

*** ADDITIONAL PARAMETERS (1 of 3) *** 

Main Rotor 

1. flapping mom of inertia 2. hub height above waterline 
3. hub fuselage station 4. hub posn right of buttline 

5. mast incidence 

Tail Rotor (enter zeros (0) if using NOTAR) 

6. height above waterline 7. hub fuselage station 
8. posn rignt of buttline 9. number of blades 

10. blade chord 11. blade radius 

12. lift curve slope 13. rotational velocity 

14. flap mom of inertia 15. delta-3 angle 

16. blade twist 

Vertical Fin 

17. height above waterline 
19. posn right of buttline 
21. CL max 

23. lift curve slope' 

0. NO CHANGES 

Input the parameter to change: 

Figure 3.2. Additional Parameters Menu, Screen 1. 


*** INPUT DATA CONTINUED (screen 2 of 8) *** 
filename 

Main Rotor 

advance ratio = mu 
inflow parameter wrt TPP « lamp 

Tip path angle « altpp*57.3 degs 
Rotor coning angle ■ ao*57.3 degs 
1st lat cyclic term-Al « Al*57.3 degs 
1st long cyclic term-Bl * Bl*57.3 degs 
lateral flapping * bls*57.3 degs 
longitudinal flapping * als*57.3 degs 
Lock number ■= lockno 


Figure 3.3. Main Rotor Calculated Output. 


All derivatives can be called from the workspace. This is described in the Output Data Instructions 
screen 1 (Fig. 3.3). The file VARLIST.TXT is found in Appendix D. 


18. fuselage station 

20. alpha zero lift 

22. dynamic pressure ratio 
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*** KEY CONTROL PARAMETERS *** 


cross coupling * xcoupls 

Designed damping 
pitch - desdmdg ft-lbs/(rad/sec) 
roll - desdrdp ft-lbs/(rad/sec) 
yaw > desdndr ft-lbs/(rad/sec) 

Control Power 
pitch b cppitch ft-lbs/in 
roll b cproll ft-lbs/in 
yaw b cpyaw ft-lbs/in 

Cooper Harper Pilot Ratings 
danqping/moment of inertia 

pitch (dM/dq) /Iyy ■ prpitch [ft-lbs/ (rad/sec) ]/ (slug ft A 2) 
roll (dR/dp)/Xxx b prroll [ft-lbs/(rad/sec)]/(slug ft A 2) 
yaw (dN/dr)/Xzz b pryaw [ft-lbs/(rad/sec)]/(slug ft A 2) 

control power/moment of inertia 
pitch (dM/in)/Iyy b cpipitch (ft-lbs/in)/(slug ft A 2) 
roll (dR/in)/Ixx • cpiroll (ft-lbs/in)/(slug ft A 2) 
yaw (dN/in)/Izz b cpiyaw (ft-lbs/in)/(slug ft A 2) 


Figure 3.4. Key Control Parameters Output Screen. 


*** OUTPUT DATA INSTRUCTIONS (screen 1 of 3) *** 

Because this subroutine generates a large number of single 
value data not shown on the output screen, a text file 
VARLIST.TXT is on this disk which lists the variable names 
for all the stability derivatives. Stability derivative 
contributions for all major aircraft components can be found 
by reading the text file VARLIST.TXT, then asking MATLAB the 
variable name corresponding to the derivative. 

Press any key to continue 


Figure 3.5. Output Data Instructions, Screen 1. 


B. TAIL ROTOR CONTRIBUTIONS 

Before stability contributions of the tail rotor can be determined, the tail rotor trim solution 
must be derived. Since only the main rotor trim solution is computed in the performance section of 
JANRAD, the tail rotor trim solution is solved in the following manner. 
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By assuming that the shaft axis of the tail rotor is perpendicular to the flight path of the 
helicopter, the tip path plane angle (a^f) is equal to longitudinal flapping (a„ x ). Additionally, no 
longitudinal or lateral cyclic pitch is introduced by the control system, hence B 1T =0 and A 1T -0. This 
reduces the inflow parameter with respect to the tip path plane from [Ref. 2:p. 189] 

- p [a„„ - (fij ♦ a,)] -A 
to 


X --- * -Cj/o — 

QR * 2p 

Tail rotor thrust T t is computed by using the main rotor torque from the main rotor trim 
solution, length of the moment arm to the tail rotor, lift of the vertical tail, and length of the moment 
arm to the vertical tail. The moment contributions of the tail rotor and vertical tail must equal the 
torque of the main rotor. This yields 


Tt = 


Cq h P A u (Q yRy ) 2 - L v l v 


lr 


Tail rotor coning is arrived at by using centrifugal forces and tail rotor thrust. This gives the 
following, where y is the Lock number of the tail rotor blades [Ref. 2:p 171]. 


a 0 = - Y 


cyo 




3 ' a (QR) 2 

Longitudinal flapping is determined by writing the flapping equations for the tail rotor in terms 
of a„, (y«r, n, and the tail rotor pitch-flap coupling angle 6 y Furthermore, this allows the lateral 
flapping to be determined in terms of the above values [Ref. 2:p 189]. 
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At a hover, there is no longitudinal or lateral flapping in the tail rotor. Because of this, the 
tail rotor equations can be resolved in terms a„, 0^, 0 IT , ht. It, and the main rotor trim conditions. 
This results in the following set of equations for the tail rotor at a hover. 



where c 4 is approximated from 



16 






Once a trim solution has been determined, stability derivatives are found in the same manner 
as for the main rotor, by perturbing about a nominal point. At a hover, the stability derivatives are 
computed in terms of the main rotor stability derivatives as in Ref. 2. 

These equations may produce flapping magnitudes that are lower than more exact solutions, 
but since the trends of the perturbations are the same in both cases, we can use these trim equations 
as the basis for solving the basic derivatives of the tail rotor [Ref. 2:p 189]. 

In addition to the tail rotor parameters input from the JANRAD opening menu, additional 
parameters are required for the stability and control analysis. Tail rotor inputs are entered in the 
Stability and Control Routine New File menu (Fig. 3.6) and changed in the Stability and Control 
Menu Additional Parameters menu, screen 1 (Fig. 3.7). 


Tail rotor height above reference datum/waterline (ft): 
Tail rotor fuselage station (ft): 

Tail rotor position right of buttline (ft): 

Number of tail rotor blades: 

Blade chord (ft): 

Tail rotor blade radius (ft): 

Average lift curve slope of tail rotor: 

Rotational velocity of tail rotor (rad/sec): 

Blade flapping moment of inertia (slug ft A 2) : 

Delta-3 angle (deg): 

Blade twist uleg) : 


Figure 3.6. Tail Rotor New File Menu. 


The program tests to see if an anti-torque device is installed by checking the number of tail 
rotor blades and the NOTAR" fan diameter. If both are zero, an error message is displayed (Fig. 
3.8). 
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*** STABILITY AMD CONTROL MENU *** 

*** ADDITIONAL PARAMETERS (1 of 3) *** 

Main Rotor 

1. flapping mom of inertia 2. hub height above waterline 
3. hub fuselage station 4. hub posn right of buttline 

5. mast incidence 

Tail Rotor (enter zeros (0) if using NOTAR) 

6. height above waterline 7. hub fuselage station 
8. posn right of buttline 9. number of blades 

10. blade chord 11. blade radius 

12. lift curve slope 13. rotational velocity 

14. flap mom of inertia IS. delta-3 angle 

16. blade twist 

Vertical Fin 

17. height above waterline 
19. posn right of buttline 
21. CL max 

23. lift curve slope' 

0. NO CHANGES 

Input the parameter to change: 

Figure 3.7. Additional Parameters Menu, Screen 1. 


You must have a tail rotor or NOTAR/thruster! 
Figure 3.8. Missing Component Warning. 


Tail rotor output is displayed in Input Variable output screen 3 and in the Tail Rotor Calculated 
Data screen (Fig 3.9). 

Tail rotor derivatives are solved for in the same manner as for the main rotor. A detailed 


listing of tail rotor derivatives is in Appendix C. 


C. FUSELAGE CONTRIBUTIONS 

The fuselage contributions to stability and control are typically arrived at through empirical 
methods, so are best estimated initially by referring to characteristics of similar helicopters. 


18. fuselage station 

20. alpha zero lift 

22. dynamic pressure ratio 
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*** INPOT DATA CONTINUED (screen 3 of 8) *** 
filename 

Tail rotor (zero if NOTAR) 

Number of blades * bt 

Blade chord * cot ft 
Blade radius * Rt ft 
Lift curve slope ■ at 
Rotational velocity * ohmt rad/sec 
Flapping moment of inertia ■ Ibt slug ft*2 
Delta-3 angle = delta3 deg 
Blade twist » thetalt deg 
Hub height above waterline * htd ft 
Hub fuselage station ■ ltd ft 
Hub position rt of buttline * ytd ft 

*** CALCULATED DATA *** 

Tail Rotor (zero if NOTAR) 

tail rotor thrust = Tt lbs 
advance ratio > mut 
inflow parameter = lampt 
Rotor coning angle * aot*57.3 degs 
lateral flapping ~ blst*57.3 degs 
longitudinal flapping = alst*57.3 degs 
Lock number * locknot 


Figure 3.9. Tail Rotor Output. 


Lift and drag moments of the fuselage depend on the angle of attack a F , and the sideslip angle 
0 r of the fuselage. Because of the downwash induced by the main rotor, a F is estimated by 


a f = 0 - y e - 


where 0 is the fuselage pitch attitude, y c is the climb angle, and 


c 




v, AqA u 


The value Vp/v,, the fuselage downwash ratio, is the only fuselage stability parameter required 
to be input into the stability and control program. It can be estimated by methods outlined in Ref. 
2. A typical value of 1.5 is recommended otherwise. 
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The other parameters input with fuselage data are the eg location, and moments of inertia of 


the aircraft. 

Input for the fuselage is made in the Stability and Control Routine New File menu (Fig. 3.10) 
and changed in the Stability and Control Menu Additional Parameters menu, screen 3 (Fig. 3.11). 


CG height above reference datum/waterline (ft): 
CG Fuselage station (ft): 

CG position right of buttline (ft): 

Ixx (slug ft A 2): 
lyy (slug ft A 2): 

Izz (slug ft A 2): 

Ixz (slug ft A 2): 

Oownwash ratio for fuselage (page 513 Prouty): 


Figure 3.10. Fuselage New File Menu. 


*** STABILITY AND CONTROL MENU *** 

*** ADDITIONAL PARAMETERS (3 of 3) *** 

CG location and Inertias/fuselage parameters 
1. eg ht. above waterline 2. eg fuselage station 
3. eg posn rt of buttline 4. Ixx 
5. Iyy 6. Izz 

7. Ixz 8. fuselage downwash ratio 

NOTAR if available (enter zeros if using tail rotor) 

9. height above waterline 10. boom fuselage station 
11. boom position left ref 12. NOTAR diameter 
13. swirl angle at boom 14. NOTAR max force 

15. thruster fuselage station 

Rigging 

16. B1 main/in defl (del e) 17. A1 main/in defl (dela) 

18. thetaOm/in defl (del c) 19. thetaOt/pedal defl (del r or p) 
20. NOTAR sleeve twist/defl 21. max rudder defl 

0. NO CHANGES 

Input the parameter to change: 


Figure 3.11. Additional Parameters Menu, Screen 3. 


Stability derivatives from the fuselage are combined with the stability derivatives of other 


components for display in the Key Control Parameters output screen. The individual contributions 
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of the fuselage to stability can be pulled from the workspace by finding the variable name in 
VARLIST.TXT and calling the value as described in the Main Rotor section. 

The values of the basic fuselage stability contributions are set in the program to typical values 
found in Ref. 2. If one desires to change these values, edit the subroutine CBODYGRP.M lines 17 
through 22, substituting desired values for the defaults. 

D. HORIZONTAL STABILIZER CONTRIBUTIONS 

Contributions of the horizontal tail are computed as for a fixed wing aircraft with one 
exception. In a helicopter, the downwash created by the rotor is very large, and in nearly all modes 
of flight, the horizontal tail is in the downwash and is affected. 

The effects of this are accounted for by adjusting the dynamic pressure on the horizontal tail 
by determining the dynamic pressure ratio q H /q, and downwash angles due to the main rotor c MH , and 
fuselage € m . These effects are generally evaluated in some empirical fashion, but some general 
guidelines follow. 

The dynamic pressure at the horizontal and vertical tail is generally lower than free stream due 
to momentum loss from the fuselage and main rotors. A first guess for q^/q can be made by using 
Figure 8.9 Ref. 2. If no information is available, the author recommends a typical value of 0.6 be 
used. 

The main rotor downwash ration for the horizontal tail has the form 

e = = v « DL 

U ’ Vj V v, 4 q 

The ratio v H /v, has a maximum theoretical value of two because of the maximum theoretical 
downwash velocity achieved by the main rotor using momentum theory, but empirical measurements 
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indicate that higher ratios are obtainable [Ref. 2]. The term v H /v, can be estimated by using Fig 8.11 
Ref. 2. The author recommends a typical value of 1.5 in the absence of more precise information. 
The fuselage down wash ratio has the form 


tm = e. 


■r*i 


de. 


\ da rj 


The fuselage down wash ratio can be estimated using Fig 8.14 Ref. 2. The fuselage down wash ratio 
includes the effects of any wings since most wings installed on a helicopter are small and do not 
warrant a separate term for downwash effects. A typical value of 0.25 is recommended when no 
other information is available. Detailed discussion on determination of these effects are contained in 
Ref. 2 page 489. 

The final effect to account for is the time lag for change in rotor downwash to reach the 
horizontal tail. This effect is accounted for through the time lag 



yielding the following equation for the angle of attack for the horizontal tail 


a 


h = 



These effects, when combined with typical computations for horizontal stabilizer stability analysis, 
yield the equations found in Appendix C. 

Input for the horizontal tail is made in the JANRAD Performance input screen with additional 
input required for stability analysis entered in Stability and Control Routine New File menu (Fig. 
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Height above reference datum/waterline (ft) 
Fuselage station (ft): 

Position right of buttline: 

Zero lift angle for horixontal tail (deg): 
Angle of incidence of horixontal tail (deg) 
Lift curve slope of horixontal tail: 
Dynamic pressure ratio (pg 489 Prouty): 
Rotor downwash ratio (pg 489 Prouty): 
Fuselage downwash ratio (pg 489 Prouty): 


Figure 3.12. Horizontal Tail New File Menu. 

3.12) and changed in the Stability and Control Menu Additional Parameters menu, screen 2 (Fig. 

3.13) . 


*** STABILITY AND CONTROL MEND *** 
*** ADDITIONAL PARAMETERS (2 of 3) *** 


Horizontal Tail 
1. height above waterline 2. 
posn right of buttline 4. 
angle of incidence €. 

dynamic pressure ratio 8 . 
fuselage downwash ratio 


Wing 

10. height above waterline 11. 
12. posn right of buttline 13. 
14. angle of incidence 15. 

16. tip cord 17, 

18. rotor downwash ratio 19. 

0. NO CHANGES 

Input the parameter to change: 


fuselage station 
alpha zero lift 
lift curve slope 
rotor downwash ratio 


fuselage station 
alpha zero lift 
lift curve slope 
root cord 

fuselage downwash ratio 


Figure 3.13. Additional Parameters Menu, Screen 2. 


Stability derivatives from the horizontal tail are combined with the stability derivatives of all 
other components for display in the Key Control Parameters output screen. The individual horizontal 
tail stability derivatives can be pulled from the workspace by finding the variable name in 
VARLIST.TXT and calling the value from the workspace. 
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E. VERTICAL STABILIZER CONTRIBUTIONS 

The stability derivatives are evaluated in the same manner as the horizontal tail. The only 
differences are that the swirl from the main rotor downwash, the induced flow from the tail rotor, 
induced velocities from the fuselage and aircraft sideslip all create an effective sideslip (Bidewash) 
angle on the vertical fin. 

If the tail rotor is mounted close to the vertical fin, the form for the downwash effect is the 
same as for the main rotor downwash on the fuselage 


tl r = 


D.L. } 



v 


where 


D.L. 




The parameter qy/q can be estimated using Figure 8.9 in Ref. 2. A typical value of 0.6 is 
recommended otherwise. 

The effect of aircraft sideslip on the vertical fin is negligible when the tail rotor is in close 
proximity to the vertical fin because the induced velocities of the tail rotor are much higher than the 
sideslip velocities typical in the linear region of this analysis. In the case where the tail rotor is far 
away from the vertical fin, or when a NOTAR” is used, a typical value of sideslip effect is assigned 
by the program in the form 
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<*»y 

df> 


dt. 

—1 - 0.06 


The induced side wash from the main rotor and other factors affecting the tail vertical fin are 
evaluated by Stability and Control section of JAN RAD and do not need to be input by the user. A 
detailed explanation of these effects in available in Ref. 2. 

Input for the vertical tail is made in the JAN RAD Performance input screen with additional 
input required for stability analysis entered in Stability and Control Routine New Hie menu (Fig. 
3.14) 


Height above reference datum/waterline (ft): 
Fuselage station (ft): 

Position right of buttline (ft): 

Zero lift angle for vertical tail (deg): 
Maximum Cl for vertical tail: 

Dynamic pressure ratio (pg 489 Prouty): 

Lift curve slope of vertical tail: 


Figure 3.14. Vertical Tail New File Menu. 


and changed in the Stability and Control Menu Additional Parameters menu, screen 1 (Fig. 3. IS). 

As with the horizontal tail, results for the vertical tail are combined with other stability 
derivatives for final output, but can be called from the workspace individually as outlined previously. 


F. WING CONTRIBUTIONS 

Most helicopters do not have wings, and the ones that do usually have wings only in the sense 
that the wings are aerodynamic hard points for weapon mounting. But compound helicopters do have 
wings that are designed to provide a large amount of lift for the helicopter. 

Stability and control analysis for a wing is the same as for a fixed wing aircraft, except that 
large amounts of rotor down wash must be accounted for. Rotor wash is accounted for in the same 
fashion as for the horizontal tail and the fuselage. The dynamic pressure at the wing is adjusted by 
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*** STABILITY AND CONTROL MENU *** 

•** ADDITIONAL PARAMETERS (1 of 3) *•* 

Main Rotor 

1. flapping mom of inertia 2. hub height above waterline 
3. hub fuselage station 4. hub poan right of buttline 

5. mast incidence 

Tail Rotor (enter zeros (0) if using NOTAR) 

6. height above waterline 7. hub fuselage station 
8. posn right of buttline 9. number of blades 

10. blade chord 11. blade radius 

12. lift curve slope 13. rotational velocity 

14. flap mom of inertia 15. delta-3 angle 

16. blade twist 

Vertical Fin 

17. height above waterline 
19. posn right of buttline 
21. CL max 

23. lift curve slope' 

0. NO CHANGES 

Input the parameter to change: 

Figure 3.15. Additional Parameters Menu, Screen 1. 


using a dynamic pressure ratio, q w /q,. The angle of attack is modified by down wash terms. 


18. fuselage station 

20. alpha zero lift 

22. dynamic pressure ratio 


e 



v hD.L. 
v, 4 q 



“ 6 * *V * a o. - Y c " t Uw - t Fn 

Input for the wing is made in the JANRAD Performance input screen with additional input 
required for stability analysis entered in Stability and Control Routine New File menu (Fig. 3.16) and 
changed in the Stability and Control Menu Additional Parameters menu, screen 2 (Fig. 3.17). 
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Height above reference datum/waterline (ft): 
Fuselage station (ft): 

Position right of buttline (ft): 

Zero lift angle for wing (deg): 

Angle of incidence of wing (deg): 

Lift curve slope of wing: 

Tip chord (ft): 

Root chord (ft): 

Rotor downwash ratio (pg 489 Prouty): 
Fuselage downwash ratio (pg 489 Prouty): 


Figure 3.16. Wing New File Menu. 


*** STABILITY AND CONTROL MEND *** 

*** ADDITIONAL PARAMETERS (2 of 3) *** 

Horizontal Tail 

1. height above waterline 

2. 

fuselage station 

3. posn right of buttline 

4 . 

alpha zero lift 

5. angle of incidence 

6. 

lift curve slope 

7. dynamic pressure ratio 

9. fuselage downwash ratio 

8. 

rotor do%mwash ratio 

Wing 

10. height above waterline 

11. 

fuselage station 

12. posn right of buttline 

13. 

alpha zero lift 

14. angle of incidence 

15. 

lift curve slope 

16. tip cord 

17. 

root cord 

18. rotor downwash ratio 

19. 

fuselage downwash ratio 

0. NO CHANGES 



Input the parameter to change: 



Figure 3.17. Additional Parameters Menu, Screen 2. 


As with the horizontal tail, results are combined with other stability derivatives for final output. 


but can be called from the workspace individually as outlined previously. 
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G. NOT AIT 

NOTAR" derivatives were obtained from proprietary information and cannot be developed 
here. The inputs for NOTAR" are entered in the Stability and Control Routine New File menu (Fig. 

3.18) and can be changed in the Stability and Control Additional Parameters Menu screen 3 (Fig. 

3.19) . NOTAR" input is echoed in the 'Input Variable’ output screen 4, but derivatives must be 
called from the workspace as described in prior sections. 


Height above reference datum/waterline (ft): 
Fuselage station (ft): 

Position right of buttline (ft): 

NOTAR boom diameter (ft): 

Swirl angle at boom(deg): 

Maximum thruster force of NOTAR (lbs): 
Thruster fuselage station (ft): 


Figure 3.18. NOTAR New File Menu. 


*** STABILITY AND CONTROL MENU *** 

*** ADDITIONAL PARAMETERS (3 of 3) *** 

CG location and Inertias/fuselage parameters 
1. eg ht. above waterline 2. eg fuselage station 
3. eg posn rt of buttline 4. Ixx 
5. Iyy 6. Izz 

7. Ixz 8. fuselage downwash ratio 

NOTAR if available (enter zeros if using tail rotor) 

9. height above waterline 10. boom fuselage station 
11. boom position left ref 12. NOTAR diameter 
13. swirl angle at boom 14. NOTAR max force 

15. thruster fuselage station 

Rigging 

16. Bl main/in defl (del e) 17. A1 main/in defl (dela) 

18. thetaOm/in defl (del c) 19. thetaOt/pedal defl (del r or p) 
20. NOTAR sleeve twist/defl 21. max rudder defl 

0. NO CHANGES 

Input the parameter to change: 


Figure 3.19. Additional Parameters Menu, Screen 3. 
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H. RIGGING 


The control input matrix [B] is automatically computed by the program, but the user must 
establish control rigging slopes and control limits. User inputs are entered in the Stability and Control 
Routine New File menu (Fig. 3.20) and can be changed in the Stability and Control Additional 
Parameters Menu screen 3 (Fig. 3.18). 

Rigging inputs are echoed in the output listings, but the effects of the rigging are best seen in 
the Key Control Parameters screen, screen view selection 3, and State Matrices Representations, 
screen view selection 3. 


Long cyclic pitch per inch defl (deg/in): 

Lateral cyclic pitch per inch defl (deg/in): 

Collective pitch per inch defl (deg/in): 

Tail rotor pitch change per inch defl or percentage of twist 
Enter 0 (zero) if using NOTAR (deg/in or deg/deg of twist): 

Maximum deflection of anti-torque from neutral for NOTAR, 
Enter 0 (zero) if using tail rotor (deg or inch travel): 

Displacement of anti-torque control until full rudder 
deflection. Enter 0 (zero) if rudder is fixed 
(deg or inch travel): 


Figure 3.20. Rigging New File Menu. 
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IV. SOFTWARE USE WITH EXAMPLES 


A. THE PROUTY HELICOPTER 

The program is initiated by a program call from the JAN RAD program. Alter entering 
MATLAB, the command JANRAD is entered and the opening menu appears (Fig. 4.1). 


Do you want to edit an existing file or create a new one? 

1. edit existing file 2. create new file >> 

Figure 4.1. Opening Menu. 

If entering a new file, JANRAD will ask for the required input. When basic data entry is 
complete, a screen appears requesting a filename for the data. Enter a six letter filename. This 
allows JANRAD to concatenate additional letters to the filename for specific output files (Fig. 4.2). 

If editing an old file, loading instructions appear (Fig. 4.3). The filename is entered as 
instructed. Data entry is performed by selecting the variable to change, then making the appropriate 


*** SAVE INSTRUCTIONS *** 

A. Save the new data to a specified file name. 

B. Do not use an extension or quotations. 

C. Use letter/number combinations of 6 characters or less. 

D. The file will be saved with a ".mat” extension. 

ex: dsgn_2 

E. If you made no changes, press < Enter >, the file will 
be saved with the original name. 

save file as: 

Figure 4.2. Save Screen. 


31 







entry. If < ENTER> is pressed without entering a value, the last number assigned to the variable 
is retained. This value is displayed on the screen during the procedure. 


*** LOAD INSTRUCTIONS *** 

A. Input the name of the file to edit. 

B. The file was saved in your previous session 
with a ".mat 1 * extension. 

C. Do not include the extension or quotations. 

ex: dsgnl 

name of input file: 

Figure 4.3. Load Instructions. 

After pressing < ENTER >, a menu appears asking which routine to run (Fig. 4.4). When 
Stability and Control is selected, three additional input screens appear for further data entry. When 
this is complete, another save menu appears. A new filename can be entered, or <ENTER> can 
be pressed to save data under the last filename. 


*** EXECUTION MENU *** 

1. Rotor Performance Analysis 

2. Stability and Control Analysis 

3. Rotor Dynamics Analysis 

4. Change data 

5. Quit 

Enter a 1, 2, 3, 4, or 5 >> 


Figure 4.4. Execution Menu. 


After data entry is complete, the program performs a stability and control analysis then prompts 
the user for desired output format (Fig. 4.5). If the user wants hard copies of the root loci or the 
control Bode plots, the plots must be viewed on the screen first to allow MATLAB to create the graph 
META files for printing. 
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Do you want the results displayed cm screen? 

NOTE: if you want a hard copy of the plots, you must 
select (1) and view them on the screen first. 

1. yes 2. no >> 

Figure 4.5. Data Display. 

Output from the program is in six basic categories as shown in Fig. 4.6. To demonstrate the 
output features, the Prouty sample helicopter at a hover [Ref. 2] will be used as an example. 


*** STABILITY AND CONTROL PROGRAM *** 

*** SCREEN VIEW MENU *** 

What do you want to see? 

1. Input data. 

2. Calculated data. 

3. State Matrices. 

4. Eigenvalues of the plants and plots of the roots. 

5. Key control parameters. 

6. Open loop transfer plots. 

0. Exit screen view. 

Enter a number; 

Figure 4.6. Screen View Menu. 

Selection 1, Input Data, is self explanatory, and will be skipped in the interest of brevity. 

Selection 2, Calculated Data, has two screens. The first screen (Fig. 4.7) gives the calculated 
main rotor parameters of interest. Note that since the helicopter is hovering, many of the quantities 
are equal to or near zero. The second screen has calculated data for the anti torque system (Fig. 4.8). 
Again note that many parameters are zero since the aircraft has no forward velocity. 

Selection 3 provides the state matrix representations of the longitudinal, lateral and coupled 
plants as well as the input matrices (Figs. 4.9, 4.10, and 4.11). The states are listed under each 
heading for convenience. These matrices can be used after exiting the program to simulate aircraft 
response to impulses, steps or any other variety of input. These matrices can also be used to develop 
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*** CALCULATED DATA (screen 

1 of 

2) *** 

prouth 



Main Rotor 



advance ratio * 

0.0 


inflow parameter wrt TPP * 

0.000 


Tip path angle » 

0.0 

degs 

Rotor coning angle * 

7 .4 

degs 

1st lat cyclic term-Al - 

0.0 

degs 

1st long cyclic term-Bl - 

-0.0 

degs 

lateral flapping « 

0.00 

degs 

longitudinal flapping > 

0.00 

degs 

Lock number * 

7.2 


press any key to continue... 




Figure 4.7. Calculated Data, Screen 1. 


*** CALCULATED DATA (screen 2 of 2)*** 
prouth 

Tail Rotor (zero if NOTAR) 


tail rotor thrust = 
advance ratio = 
inflow parameter = 
Rotor coning angle = 
lateral flapping = 
longitudinal flapping = 
Lock number * 


1349.5 lbs 
0.0 
-0.073 

-0.1 degs 
0.00 degs 
0.00 degs 
2.7 


press any key to continue... 


Figure 4.8. Calculated Data, Screen 2. 


response to impulses, steps or any other variety of input. These matrices can also be used to develop 
state feedback control for the aircraft. 


Selection 4 displays the eigenvalues and the root loci of the uncoupled and coupled plants. 
Before doing this, an instruction screen is displayed to assist in understanding the format of the output 
(Fig. 4.12). 

The second screen displays the longitudinal and the lateral uncoupled eigenvalues, then 
immediately plots the roots (Fig. 4.13). The plots are META filed by MATLAB after display. To 
recover them on a hard copy, one must use the Graphics Post Processing (GPP) utility after exiting 
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*** CALCULATED DATA (screen 

1 of 

2) *** 

prouth 



Main Rotor 



advance ratio ■ 

0.0 


inflow parameter wrt TPP ■ 

0.000 


Tip path angle * 

0.0 

degs 

Rotor coning angle * 

7.4 

degs 

1st lat cyclic term-Al * 

0.0 

degs 

1st long cyclic term-Bl m 

-0.0 

degs 

lateral flapping = 

0.00 

degs 

longitudinal flapping * 

0.00 

degs 

Lock number « 

7.2 


press any key to continue... 





Figure 4.7. Calculated Data, Screen 1. 


*** CALCULATED DATA (screen 2 of 2)*** 
prouth 

Tail Rotor (zero if NOTAR) 

tail rotor thrust * 1349.S lbs 
advance ratio = 0.0 

inflow parameter = -0.073 
Rotor coning angle = -0.1 degs 

lateral flapping = 0.00 degs 

longitudinal flapping = 0.00 degs 

Lock number « 2.7 

press any key to continue... 


Figure 4.8. Calculated Data, Screen 2. 

response to impulses, steps or any other variety of input. These matrices can also be used to develop 
state feedback control for the aircraft. 

Selection 4 displays the eigenvalues and the root loci of the uncoupled and coupled plants. 
Before doing this, an instruction screen is displayed to assist in understanding the format of the output 
(Fig. 4.12). 

The second screen displays the longitudinal and the lateral uncoupled eigenvalues, then 
immediately plots the roots (Fig. 4.13). The plots are META filed by MATLAB after display. To 
recover them on a hard copy, one must use the Graphics Post Processing (GPP) utility after exiting 
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Longitudinal uncoupled plant (A or F depending on notation) 
States are [u w q theta] 


-0.0106 

0 

2.5068 

-32.2000 

0 

-0.2698 

0 

0.0002 

0.0038 

0 

-0.8878 

0 

0 

0 

0.9977 

0 


Longitudinal uncoupled input matrix (B or G depending on 
notation) 

Inputs are [longitudinal cyclic, collective, lateral cyclic. 


pedals] 




1.2248 

- 0.0000 

-0.0738 

0 

0 

-5.8773 

0 

0 

-0.4338 

0.0000 

0.0261 

•0.0238 

0 

0 

0 

0 


press any key to continue... 


Figure 4.9. Longitudinal State Matrix Representation. 


Lateral/directional uncoupled plant (A or F 

depending on 

notation) 

States are 

[v p phi r psi] 



-0.0305 

-2.6260 32.1267 

0.7353 

0 

-0.0064 

-1.0274 0 

0.0783 

0 

0 

1.0000 0 

• 0.0000 

0 

0.0130 

0.0783 0 

-3.3360 

0 

0 

0 0 

1.0000 

0 

Lateral/directional uncoupled input matrix (B or G depending 

on notation) 



Inputs are 
pedals] 

[longitudinal cyclic. 

collective. 

lateral cyclic. 

0.1200 

0.0000 0.7533 

-1.9736 


0.0486 

0.0000 0.3049 

•0.2101 


0 

0 0 

0 


0 

0.3001 0 

1.2959 


0 

0 0 

0 


press any key to continue_ 




Figure 4.10. Lateral/Directional State Matrix Representation. 


the program as explained in the closing screen. 

The first root in the lateral directional plant is an artifact of the system representation and has 
no physical meaning. The other modes are the Dutch-roll and spiral modes. The longitudinal modes 


shown are the short period and phugoid. 
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Coupled plane (A or F depending on notation) 
States are (u w q theta v p phi r psi] 


Columns 

1 through 

7 





-0.0106 

0 

2.S068 

-32.2000 

-0.0057 

-0.7469 

0 

0 

•0.2698 

0 

0.0002 

0 

0 

2.1710 

0.0038 

0 

-0.8878 

0 

0.0013 

0.2645 

0 

0 

0 

0.9977 

0 

0 

0 

0 

0.0057 

0 

-0.7469 

•0.0000 

-0.0305 

-2.6260 

32.1267 

0.0023 

0 

-0.3023 

0 

-0.0064 

-1.0274 

0 

0 

0 

0.0000 

0 

0 

1.0000 

0 

0 

-0.0198 

0 

0 

0.0130 

0.0783 

0 

0 

0 

0 

0 

0 

0 

0 


Columns 8 through 9 
0 0 

0 0 


0 0 

0 0 

0.0266 0 

0.0674 0 

0.7353 0 

0.0783 0 

- 0.0000 0 

-3.3360 0 

1.0000 0 


Coupled input matrix (B or G depending on notation) 

Inputs are [longitudinal cyclic, collective, lateral cyclic, 
pedals] 


1.2248 

-0.0000 

-0.0738 

0 

0 

-5.8773 

0 

0 

-0.4338 

0.0000 

0.0261 

-0.0238 

0 

0 

0 

0 

0.1200 

0.0000 

0.7533 

-1.9736 

0.0486 

0.0000 

0.3049 

-0.2101 

0 

0 

0 

0 

0 

0.3001 

0 

1.2959 

0 

0 

0 

0 


press any key to continue... 


Figure 4.11. Coupled State Matrix Representation. 


Comparing the root estimates of the longitudinal and lateral plants with the roots of the coupled 
plant in the third screen shows close correlation. 

The plots of the eigenvalues are in the Argand plane (Figs 4.15, 4.16, and 4.17). Note that 
the plant has some unstable roots. This is typical of helicopters. Very rarely are all the modes stable 
without some form of stability augmentation. 
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After you view the root loci plot, a meta file is made. 

When you are done a screen will tell you the file names 
of the meta files. To get a hard copy of the plots, you 
must graphics post process (GPP) the files for your 
particular printer set-up then, print. 

NOTE: If ALL roots are real, MATLAB will NOT plot them 

in the Argand plane, but will plot the root against 
its position in the vector (e.g. the first root would 
be plotted as (l.root)) 

press any key to continue... 


Figure 4,' 1. Eigenvalue and Root Loci Instruction Screen. 


*** EIGENVALUES 

+ ** 


prouth 



Uncoupled 



Longitudinal plant 



Root 

wn 

damping 

-1.0155 

1.0155 

1.0000 

0.0585 + 0.3398i 

0.3448 

-0.1698 

0.0585 - 0.3398i 

0.3448 

-0.1698 

-0.2698 

0.2698 

1.0000 

Lateral/Directional 

plant 


Root 

wn 

damping 

0 

0 

- 1.0000 

-3.3433 

3.3433 

1.0000 

-1.1784 

1.1784 

1.0000 

0.0639 + 0.4025i 

0.4075 

-0.1568 

0.0639 - 0.4025i 

0.4075 

-0.1568 

press any key to continue... 



Figure 4.13. Uncoupled Eigenvalues. 


The final screen for this section (Fig. 4.18) tells the user the filenames for the graphs. The 


user can write these down, or refer to Appendix E for the filenames when bard copies are desired. 


The Key Control Parameters screen 2 (Figs 4.19 and 4.20), selection 3, provides numbers that 


can be used to determine Cooper-Harper pilot ratings and compliance with MIL-H-8501. The screen 


also tells the amount of cross coupling between flapping modes of the main rotor system. Though 
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Coupled Plant 


*** BXGKNVALUBS *** 
prouth 


Root wn damping 


0 



0 

- 1.0000 

-3.3432 



3.3432 

1.0000 

-1.0970 

♦ 

0.2583i 

1.1270 

0.9734 

-1.0970 

- 

0.2583i 

1.1270 

0.9734 

0.0760 

♦ 

0.3983i 

0.4054 

-0.1873 

0.0760 

- 

0.3983i 

0.4054 

-0.1873 

0.0461 

♦ 

0.3582i 

0.3611 

-0.1276 

0.0461 

- 

0.3582i 

0.3611 

-0.1276 

-0.2689 



0.2689 

1.0000 


press any key to continue... 


Figure 4.14. Coupled Eigenvalues. 
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Figure 4. IS. Longitudinal Roots. 


Plots are saved under the following filenames: 

Longitudinal roots - rootlon.met 
Lateral/Directional roots - rootlat.met 
Coupled roots - rootcoup.met 

press any key to continue... 


Figure 4.18. Eigenvalues Closing Screen. 
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Figure 4.17. Coupled Roots. 


MIL-H-8501 has been superseded by ADS-33C, the numbers still provide information on the damping 
characteristics and control power of the helicopter. 

The last selection, open loop transfer plots, provides the transfer function from a control input 
to a state output. The purpose of these plots is to allow the designer to design automatic flight control 
systems and analyze the best techniques for controlling unstable modes. For example, to see the 
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*** KEY CONTROL PARAMETERS (scrssn 1 of 2) *** 
prouth 

cross coupling * 0.08 

Designed damping 

pitch - -35513.5 ft-lbs/(rad/sec) 
roll « -35957.8 ft-lbe/(rad/sec) 
yaw • -H6759.1 ft-lbs/(rad/sec) 

Control Power 

pitch ■ -17351.2 ft-lbs/in 
roll - 10672.0 ft-lbs/in 
yaw « 45355.3 ft-lbs/in 

press any key to continue... 


Figure 4.19. Key Control Parameters, Screen 1. 


*** KEY CONTROL PARAMETERS (screen 2 of 2) •** 
prouth 

Cooper Harper Pilot Ratings 
damping/moment of inertia 

pitch (dM/dq)/Iyy * -0.89 [ft-lbs/(rad/sec)]/(slug ft*2) 

roll (dR/dp)/Ixx * -1.03 [ft-lbs/(rad/sec)1/(slug ft*2) 

yaw (dN/dr)/Izz » -3.34 [ft-lbs/(rad/sec)]/(slug ft A 2) 

control power/moment of inertia 
pitch (uM/in)/Iyy * -0.43 (ft-lbs/in)/(slug ft*2) 

roll (dR/in)/Ixx * 0.30 (ft-lbs/in)/(slug ft~2) 

yaw (dN/in)/Izz « 1.30 (ft-lbs/in)/(slug ft“2) 

press any key to continue... 


Figure 4.20. Key Control Parameters, Screen 2. 


effects of collective pitch on forward velocity, u, select the longitudinal set of plots from the screen 


(Fig. 4.21) and page through the plots until you get to "Collective Pitch to U." It can be seen as 


shown in Fig. 4.22, that there is little control gain from collective to forward velocity at a hover. 


Obviously, collective would be a poor choice for controlling fore and aft motion at a hover. 


If one looks at "Longitudinal Cyclic to U," it can be seen that there is high gain from input 


to output, especially at the short period frequency (Fig. 4.23). This shows that longitudinal cyclic 


is probably a good choice for controlling fore and aft motion. 
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After you view a bode plot of the transfer function from 
input to state output, a meta file is made. When you exit, 
a screen will tell you the file names of the meta files. 

To get a hard copy of the plots, you must graphics post 
process (GPP) the files for your particular printer set-up 
then print. 

While viewing a plot, press any hey to go to the next plot 
Do you want to see longitudinal or lateral/directional plots? 

1. Longitudinal (eight plots total). 

2. Lateral Directional (ten plots total). 

Enter a number : 


Figure 4.21. Open Loop Transfer Plots, Screen 1. 



Figure 4.22. Collective to Forward Velocity. 


Another observation from this plot is that at low frequencies (e.g. when you hold the cyclic 
steady) there is good control response. This is consistent with flight experience. When the cyclic 
is displaced forward, the aircraft pitches forward as the cyclic is moved, until some equilibrium pitch 
attitude is reached. When the cyclic is stopped, the attitude remains fixed and the aircraft moves 
forward. Likewise, when the cyclic is moved back and forth very rapidly, the aircraft may vibrate, 
but the pitch attitude remains essentially fixed because the aircraft cannot move as quickly as the 










Figure 4.23. Longitudinal Cyclic to Forward Velocity. 


control input, hence, the aircraft does not change its average velocity. 

These effects can be seen in the control responses of all the controls. In the lateral-directional 
plant, the response from lateral cyclic to sideslip, v, has the same type of peak near the Dutch-roll 
mode, but response above the Dutch-roll frequency is reduced (Fig. 4.24). 



Figure 4.24. Lateral Cyclic to Lateral Velocity (Sideslip). 
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A complete set of Bode plots for the prouty helicopter can be found in Appendix G. The 
shapes of these curves is typical of helicopters with single main rotors and a tail rotor. 

Hard copies of the plots are obtained in the same manner as the root locus plots. A complete 
list of command response META files is located in Appendix E. 

The output from the screens is written to a diary after exiting the Screen View menu. The 
Output Data Instruction menus in Figs. 4.25 through 4.27 explain the data recovery procedures. 


*** OUTPUT DATA INSTRUCTIONS (screen 1 of 3) *** 

Because this subroutine generates a large number of single 
value data not shown on the output screen, a text file 
VARLIST.TXT is on this disk which lists the variable names 
for all the stability derivatives. Stability derivative 
contributions for all major aircraft components can be found 
by reading the text file VARLIST.TXT, then asking MATLAB the 
variable name corresponding to the derivative. 


Figure 4.25. Output Data Instructions, Screen 1. 


*** OUTPUT DATA INSTRUCTIONS (Screen 2 of 3) *** 

A. Data from the output screen saved to a file named: 

"fi1ename1.stb" 

B. This is a text file, use the TYPE command to view the file 
or use a text editor to view/print the file. 

C. Matrix and vector data saved to a default file named: 
"mstabdat.mat" 

D. This is a ".mat" binary file, use the LOAD command to 
retrieve the data for plotting. 

E. Rename "mstabdat.mat" to another ".mat" file. 

The file "mstabdat.mat" will be overwritten when 
the program is executed. 

F. Do not rename the file as "filenamel.mat" 

The file "filenamel.mat" is already on disk 
and used for future editing. 


Figure 4.26. Output Data Instructions, Screen 2. 
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*** OUTPUT DATA INSTRUCTIONS (Bcraen 3 of 3) *** 

A. Single value data saved to a default file named: 
"vstabdat.mat” 

B. This ia a ".mat" binary file, use the LOAD command to 
retrieve the data for plotting. 

C. Rename "vstabdat.mat" to another ".mat* file. 

The file "vstabdat.mat" will be overwritten when 
the program is executed. 

D. Do not rename the file as "filenamel.mat" 

The file "filenamel.mat" is already on disk 
and used for future editing. 

*** END STABILITY AND CONTROL ROUTINE *** 


Figure 4.27. Output Data Instructions, Screen 3. 

B. A HELICOPTER WITH NOTAR" 

A final example will be a notional helicopter designed for the 1993 annual American Helicopter 
Society (AHS) Design Competition, the Arapaho. The forward flight mode of the helicopter will be 
used. Since the procedures for entering and displaying data have been covered, only the salient points 
of the design will be addressed. 

The helicopter incorporated a NOTAR" tailboom. As a result, the tail rotor derivatives are 
unusual as are the command transfer Bode plots. 

The Output file will have a zero displayed for all tail rotor computations except the tail rotor 
thrust. The thrust displayed in this case will be the NOTAR" thruster force (Fig 4.28). 

By comparing the longitudinal control matrices with and without a tail rotor (Fig. 4.10 and 
4.29) it can be seen that a change in tail rotor thrust has a pitching moment contribution (look at the 
B(3,4) term). Because the tail rotor rotates about an axis that is parallel to the Y control axis, this 
contribution is to be expected with a change in tail rotor thrust requirement. In contrast, the 
NOTAR" does not have this contribution because the thruster fan rotates on an axis parallel to the X 
control axis. 
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*** CALCULATED DATA (screen 2 of 2)*** 
arapac 

Tail Rotor (zero if NOTAR) 


tail rotor thrust 
advance ratio 
inflow parameter 
Rotor coning angle 
lateral flapping 
longitudinal flapping 
Lock number 


264.2 lbs 
0.0 
0.000 
0.0 degs 
0.00 degs 
0.00 degs 
0.0 


press any key to continue... 


Figure 4.28. Tail Rotor Output. 


Longitudinal uncoupled plant (A or F depending on notation) 

States are [u w q theta] 

-0.0363 

0.0567 3.3538 

-32.1459 

0.0685 

0.3972 84.6890 

-1.8657 

-0.0098 

-0.0204 -0.5489 

0 

0 

0 0.9999 

0 

Longitudinal uncoupled input matrix (B or 6 depending on notation) 
Inputs are [longitudinal cyclic, collective, lateral cyclic, pedals] 

7.0005 

-0.5855 -2.0435 

0 

26.0454 

•10.9228 0 

0 

-1.2115 

0.5179 0.3537 

0 

0 

0 0 

0 

press any key to continue... 



Figure 4.29. Arapaho Control Matrix. 


Figure 4.30 demonstrates the difference in flight control response with airspeed. In the hover 
case, longitudinal cyclic has little direct effect on vertical velocity. But notice that in the forward 
flight regime, a longitudinal deflection has a great effect on vertical velocity. A zoom climb can be 
done in forward flight but not at a hover. 


45 







Gain (dB) 


Open loop response Longitudino! Cyclic to W, Cruise 



Figure 4.30. Longitudinal Cyclic to Vertical Velocity. 
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V. CONCLUSIONS AND RECOMMENDATIONS. 


A. CONCLUSIONS 

The aerodynamicist should not be intimidated by the complexity of helicopters. The 
complexities of determining the stability characteristics of a helicopter are formidable, but helicopter 
flight dynamics are essentially no different than fixed wing flight dynamics. Modern software reduces 
helicopter analysis to a routine procedure, making the analysis of these characteristics much less time 
consuming. 

B. RECOMMENDATIONS 

To improve the fidelity of JANRAD Stability and Control, the following actions are 
recommended. 

a) Expand the state matrix to include a,. M , b ltM , Cya M , Ch/o M , Cq/o M , fi, X’ a„ T , b llT , Cj/o T , 
Ch/o T , Cq/o T , [Lj and Xj for the tail and main rotors. This representation would allow students to run 
simulations on MATLAB that could show rotor flapping, cyclic pitch changes, and fluctuations in the 
coefficients of torque, thrust and H-force. 

b) Expand the input screens to allow selecting different lift and drag characteristics for the 
fuselage. Even though this type of information is derived empirically, studies of existing helicopter 
with appropriate flight test data could be done without the added artificiallity of a particular fuselage 
configuration. JANRAD currently uses fuselage aerodynamic characteristics of the Prouty sample 
helicopter as outlined in section III. C. 

c) Improve the fidelity of the tail rotor trim solution. The method used to derive the tail rotor 
trim solution is extremely rough compared to the trim solution of the main rotor. 
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d) Expand the wing inputs to include dihedrally mounted wings and wing sweep. Even though 
most helicopters do not use wings for lift or control, compound helicopters are likely to require a 
higher fidelity wing solution. 
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APPENDIX A. NOTATION. 


Symbol 

Definition 

Units 

a 

Lift curve slope 

1/rad 

*0 

Coning angle 

rad, deg 


First harmonic of longitudinal 
flapping wrt shaft axis 

rad, deg 

b 

Number of blades 


b„ 

First harmonic of lateral 
flapping wrt shaft axis 

rad, deg 

c 

Chord 

ft 


Coefficient of drag 


eg 

Center of gravity 


C, 

Coefficient of lift 


C m 

Coefficient of moment 


e 

Hinge offset 

ft 

f 

Equivalent flat plate area 

ft 2 

g 

Acceleration due to gravity 

ft/sec 2 


Incidence 

rad, deg 

1 

Fuselage station 

ft 

P 

Roll rate 

rad/sec 

q 

Pitch rate 

rad/sec 

q 

Dynamic pressure 

ib/ft 

r 

Blade element radius 

ft 

r 

Yaw rate 

rad/sec 

u 

forward velocity 

ft/sec 

u 

Input vector 


V 

Induced velocity 

ft/sec 

V 

Lateral velocity 

ft/sec 

w 

Vertical velocity 

ft/sec 

y 

Output vector 


y 

Buttline position 

ft 

z 

Waterline position 

ft 

A 

Plant matrix 


A 

Area of disk 

ft 2 

\ 

Blade area 

ft 2 

A, 

First lateral harmonic of 
blade feathering 

rad,deg 

B 

Input matrix 


B, 

First longitudinal harmonic of 
blade feathering 

rad,deg 

C 

Output matrix 


C h 

Coefficient of horizontal force 



Coefficient of power 
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C, 

Coefficient of torque 


C, 

Coefficient of thrust 


D 

Direct transition matrix 


D.L. 

Disk loading 

lb/ft 2 

D 

Drag 

lb 

GW 

Gross weight 

lb 

H 

Horizontal force on rotor 

lb 

U 

Flapping moment of inertia 

slug ft 2 


Rolling moment of inertia 

slug ft 2 

by 

Pitching moment of inertia 

slug ft 2 


Yawing moment of inertia 

slug ft 2 

L 

Lift 

lb 

M 

Pitching moment 

ft-lb 

N 

Yawing moment 

ft-Ib 

Q 

Torque 

ft-lb 

R 

Blade radius 

ft 

R 

Rolling moment 

ft-lb 

T 

Thrust 

lb 

X 

State vector 


a 

Angle of attack 

rad, deg 

0 

Blade flapping angle 

rad, deg 

0 

Sideslip angle 

rad, deg 

y 

Lock number 


7c 

Climb angle 

rad, deg 


Pitch-flap coupling angle 

rad, deg 

c 

Down wash angle 

rad, deg 

V 

Sidewash angle 

rad, deg 

e 

Pitch angle 

rad, deg 

e 

Blade pitch 

rad, deg 

e, 

Blade twist 

rad, deg 

X 

Inflow ratio wrt swashplate 


V 

Inflow ratio wrt tip path plane 


M 

Tip speed ratio 


p 

Density of air 

slug/ft 3 

a 

Rotor Solidity 


4> 

Roll angle 

rad, deg 

0 

Inflow angle 

rad 

¥ 

Yaw angle 

rad, deg 

0 

Rotational velocity of rotor 

rad/sec 

QR 

Rotor tip speed 

ft/sec 
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Subscript 

Definition 

0 thru 1 

Blade position in r/R 

0 

Nominal 

0 

Reference 

b 

Per blade 

i 

Induced 

F 

Fuselage 

H 

Horizontal stabilizer 

M 

Main rotor 

N 

NOTAR" 

Q 

Torque 

T 

Tail rotor 

T 

Thrust 

TPP 

Tip path plane 

V 

Vertical stabilizer 
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APPENDIX B. NASA STATE SPACE REPRESENTATION. 
NASA state space representation from Ref. 3. 















APPENDIX C. SUMMARY OF STABILITY DERIVATIVES. 


A summary of stability derivatives extracted from Ref. 2. 
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Main Rotor Derivatives near Mover 
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Basic Tall Rotor Derivatives In Forward Flight 
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Nondlmenslonel Horizontal Stabilizer Derivatives 
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Nondimensional Vertical Stabilizer Derivatives 
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Vertical Stabilizer Derivatives in Forward Flight 
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APPENDIX D. VARIABLE LIST. 

Variable list for JANRAD Stability Routine 

A - rotor disk area 

Ab - area of blades 

Abt - area of tail rotor blades 

Ap • projected area of fuselage under rotor 

A1 - lateral cyclic pitch 

als - longitudinal flapping 

a - main rotor blade lift curve slope 

Ah - area of the horizontal stab 

ah - lift curve slope of the horizontal stab 

alph - angle of attack of horizontal stab 

alphas • angle of attack WRT shaft axis 

alploh - alpha-zero-lift of the horizontal stab ( at the trim condition ) 

alplov - alpha-zero-lift of the vertical stab 

alplow - alpha-zero-lift of the wing 

alpv - angle of attack of vertical stab 

alpw - angle of attack of wing 

altpp - mean alpha TPP 

ao - main rotor coning angle 

aot - tail rotor coning angle 

At - disc area of the tail rotor 

at - tail rotor blade lift curve slope 

at - tail rotor blade lift curve slope 

Av - area of the vertical stab 

av - lift curve slope of the vertical stab 

Aw - area of the wing 

aw - lift curve slope of wing 

B1 - longitudinal cyclic pitch 

bis -lateral flapping 

beta - trim sideslip angle 

bh - span of the horizontal stab 

bv - semispan of the vertical stab 

bw - span of the wing 

c - main rotor chord 

cdoh - Cdo of the horizontal stab 

cdow - Cdo of the wing 

ch - coefficient of ???? 

chsig - Ch/sigma 

clvertmax - maximum cl of the vertical stab with full rudder defl 

emu - circulation coefficient for NOTAR" 

emun - design Cmu for NOTAR" boom in hover (.3-.6) 

cot - tail rotor chord 

cpipitch - Control power/inch pitch 


68 







cpiroll - Control power/inch pitch 
cpiyaw - Control power/inch pitch 
cppitch - Control power pitch 
cproll - Control power roll 
cpyaw - Control power yaw 
ct • coefficient of thrust 
ctsig - Ct/sigma 

ctsigt - Ct/sigma of the tail rotor 
ctr - wing root chord 
ctw - wing tip chord 
cq - coefficient of torque 
cqsig - Cq/sigma 

deldv - biplane effect (mutual interference of tail rotor and vertical stab) 

delih - span efficiency factor of the horizontal stab 

deliv - span efficiency factor of the vertical stab 

deliw - span efficiency factor of the wing 

delta3 - tail rotor delta-3 angle 

dclvmax - maximum rudder deflection on the vertical stab 

desdmdq - designed pitch damping (dM/dq) 

desdndr - designed yaw damping (dN/dr) 

desdrdp - designed roll damping (dR/dp) 

detafdalpfh - fuselage downwash ratio for the horizontal stab 

detafdalpfw - fuselage downwash ratio for the wing 

detafdbeta - fuselage sideslip ratio (= depsilondalpha) 

dian - NOTAR” boom diameter 

e - blade offset 

GW - gross weight 

g - acceleration due to gravity 

hh - horizontal stab vertical offset 

hhd - height from waterline to horizontal tail 

hm - main rotor vertical offset 

hmd - height from waterline to main rotor hub 

hslot - total of slot heights for NOTAR” boom 

ht - tail rotor vertical offset 

htd - height from waterline to tail rotor hub 

htn - NOTAR" thruster vertical offset 

hv - vertical stab vertical offset 

hvd - height from waterline to vertical fin 

hw - wing vertical offset 

hwd - height from waterline to wing 

htnd - height from waterline to NOTAR" 

lb - blade flapping inertia 

Ibt - tail rotor blade flapping moment of inertia 

Ixx - mass moment of inertia about x axis 

Iyy - mass moment of inertia about y axis 

Ixz • mass moment of inertia about xz plane 

Izz - mass moment of inertia about z axis 

Ic - Ixx*Izz 
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i - shaft incidence main rotor 

ih • angle of incidence of horizontal stab 

iw - angle of incidence of wing 

lamp - lambda’ - inflow ratio WRT TPP 

Ih - horizontal stab longitudinal offset 

Utd - fuselage station of horizontal tail 

lm - main rotor longitudinal offset 

lmd • fuselage station of main rotor hub 

It - tail rotor longitudinal offset 

ltd - fuselage station of tail rotor hub 

Itn • NOTAR" thruster longitudinal offset 

Itnd - fuselage station of NOTAR" boom 

lttnd - fuselage station of NOTAR" thruster 

Iv - vertical stab longitudinal offset 

Ivd - fuselage station of vertical tail 

iw - wing longitudinal offset 

Iwd - fuselage station of wing 

lockno - main rotor lock number 

locknot - tail rotor lock number 

It - tail rotor longitudinal offset 

Iv - vertical stab longitudinal offset 

maxr - maximum rudder deflection 

mu - advance ratio 

mut - tail rotor advance ratio 

n - a counter in CTPLOTS.M 

nt * number of tail rotor blades 

ohm - omega 

ohmt - tail rotor omega 

phin - NOTAR” thruster sleeve rotation angle 

pho - roll trim attitude (euler angle) 

prpitch - Pilot rating (dM/dq)/Iyy) 

prroll - Pilot rating (dR/dp)/Ixx) 

pryaw - Pilot rating (dN/dr)/Izz) 

q - dynamic pressure, 1/2 rho V*2 

qhq - horizontal tail dynamic press ratio (qh/q) 

qvin - dynamic press due to downwash at the NOTAR” slots 

qvq - vertical tail dynamic press ratio (qv/q) 

R - blade radius 

Rt - tail rotor blade radius 

rho - air density 

sidearm - maximum pedal travel or twist grip deflection 

sigma - solidity 

sigmat - tail rotor solidity 

swirl • mean swirl angle off main rotor measured at the NOTAR” slots 
T - thrust 

Tt - thrust of the tail rotor 
thetal - main rotor blade twist 
thetalt - tail rotor blade twist 



theUo - blade pitch at root 
thetat - induced angle at the tip 
tho - trim pitch attitude (euler angle) 

V - forward velocity 

vl - induced flow velocity at a hover 

vlf - approximation for vl at forward airspeed 

vfvl - rotor down wash ratio for the fuselage 

vhvl - rotor down wash ratio for horizontal stab 

vwvl - rotor downwash ratio for wing 

uo - x initial velocity 

wo - y initial velocity 

wo - z initial velocity 

xcg - height from waterline to eg 

xcouple • designed cross coupling 

Ytmaxn - maximum Y force from NOTAR” thruster 

yeg - length from buttline to eg 

Yv - side-force of vertical stab 

yhd - length from buttline to horizontal tail 

ym - main rotor lateral offset 

ymd - length from buttline to main rotor hub 

ytd - length from buttline to tail rotor hub 

ytnd - length from buttline to NOTAR" 

ywd - length from buttline to wing 

yvd - length from buttline to verticle tail 

zcg - fuselage station of eg 

Zh - vertical force of horizontal stab 

Zw - vertical force of wing 

CRUISE 

CRUISE -- BASIC MR Derivatives 

dalda - dals/dAl 
daldb - dals/dBl 
dal diamp - dals/d(lambda') 
daldmu - dals/dmu 
daldp - dals/dp 
daldq - dals/dq 
daldtheto - dals/dtheta(o) 
dbetdydot - d(Beta)/dYdot 
dblda - dbls/dAl 
dbldb - dbls/dBl 
dbldlamp - dbls/d(lambda’) 
dbldmu • dbls/dmu 
dbldp - dbls/dp 
dbldq - dbls/dq 
dbldtheto - dbls/dtheta(o) 
dchsigda - d(Ch/sigma)/dals 
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dchsigdb - d(Ch/sigma)/dbls 
dchsigdlamp • d(Ch/sigma)/d(lambda’) 
dchsigdmu • d(Ch/sigma)/dmu 
dchsigdtheto - d(Ch/sigma)/dtheta(o) 
dctsigdiamp - d(Ct/sigma)/d(lambda’) 
dctsigdmu - d(Ct/sigma)/dmu 
dctsigdtheto - d(Ct/sigma)/dtheta(o) 
dcqsigdlamp - d(Cq/sigma)/d(lambda’) 
dcqsigdmu - d(Cq/sigma)/dmu 
dcqsigdtheto - d(Cq/sigma)/dtheta(o) 
dcysigdb - d(Cy/sigma)/dbls 
dlampdxdot • d(lambda’)/dXdot 
dlampdydot - d(lambda’)/dYdot 
dlampdzdot - d(lambda’)/dZdot 
dmdalsm - dM/dals - main rotor stiffness 
dmudxdot - dmu/dxdot 
drdblsm - dR/dBls - main rotor stiffness 

CRUISE - MR derivatives 

dmdalm - dM/dAl 
dmdbJm - dM/dBl 
dmdpm - dM/dp 
dmdqm - dM/dq 
dmdthetom - dM/dtheta(0) 
dmdxdotm - dM/dXdot 
dmdydotm - dM/dYdot 
dmdzdotm - dM/dZdot 
dndrm - dN/dr 
dndtheom - dN/dtheta(0) 
dndxdotm - dN/dXdot 
dndzdotm - dN/dZdot 
drdalm - dR/dAl 
drdblm - dR/dBl 
drdpm - dR/dp 
drdqm - dR/dq 
drdthetom - dR/dtheta(0) 
drdxdotm * dR/dXdot 
drdydotm - dR/dYdot 
drdzdotm - dR/dZdot 
dxdalm - dX/dAl 
dxdblm - dX/dBl 
dxdpm - dX/dp 
dxdqm - dX/dq 
dxdthetom - dX/dthetaO 
dxdxdotm - dX/dXdot 
dxdydotm - dX/dYdot 
dxdzdot - dX/dZdot 
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dydalm - dY/dAl 
dydblm - dY/dBl 
dydpm • dY/dp 
dydqm - dY/dq 
dydthctom - dY/dtheta(0) 
dydxdotm - dY/dXdot 
dydydotm - dY/dYdot 
dydzdotm - dY/dZdot 
dzdblm - dZ/dbls 
dzdrm - dZ/dr 
dzdthetom - dZ/dtheta(0) 
dzdxdotm - dZ/dXdot 
dzdzdotm - dZ/dZdot 

CRUISE - TAIL ROTOR Derivatives 

dctsigdlampt - d(Ct/sigma)/d(lambda’) 

detsigdmut - d(Ct/sigma)/d(mu) 

detsigdthetot - d(Ct/sigma)/d(thetaO) 

dlampdydott - d(lambda’)/dYdot 

dydxdott - dY/dXdot 

dydydott - dY/dYdot 

dydpt - dY/dp 

dydrt - dY/dr 

dydthetot - dY/d(thetaO) 

drdxdott - dR/dXdot 

drdydott - dR/dYdot 

drdpt - dR/dp 

drdrt- dR/dr 

drdthetot - dR/dthetaO) 

dndxdott - dN/dXdot 

dndydott - dN/dYdot 

CRUISE -- Vi AL STAB DERIVATIVES 

dalpvdxdotv - d(alpha vert)/dXdot 

dalpvdydotv - d(alpha vert)/dYdot 

dbetadydotv - d(beta)/dYdot 

detatvdxdotv - d(eta tail rotor)/dXdot 

detatvdydotv - d(eta tail rotor)/dYdot 

detafdydotv - d(eta fuselage)/dYdot 

drdxdotv - dR/dXdot 

drdydotv - dR/dYdot 

dxdxdotv - dX/dXdot 

dxdydotv - dX/dYdot 

dydxdotv - dY/dXdot 

dydydotv - dY/dYdot 

dydpv - dY/dp 
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dydrv - dY/dr 

CRUISE - HORIZONTAL STAB DERIVATIVES 

dalphdxdoth - d(alpha horiz)/dXdot 
dalphdzdbldoth - d(aJpha horiz)/dZ double dot 
dalphdzdoth - d(alpha horiz)/dZdot 
detamhdxdoth - d(eta main rotor)/dXdot 
detafbdzdotb - d(eta fuse)/dZdot 
detamhdzdoth - d(eta main rotor)/dZdot 
dgamdzdoth - d(gamma)/dZdot 
dxdxdoth • dX/dXdot 
dxdzdbldoth - dX/dZ double dot 
dxdzdoth - dX/dZdot 
dzdxdotb - dZ/dXdot 
dzdzdbldoth - dZ/dZ double dot 
dzdzdoth - dZ/dZdot 

CRUISE ~ FUSELAGE 

the following are from the curves in Appendix A of Ref. 2: 

dfdalpf - d(D/q)/d(alpha fuse) 
dlqdalpf - d(L/q)/d(alpha fuse) 
dmqdalpf - d(M/q)/d(alpha fuse) 
dnqdbetaf - d(N/q)/d(beta) 
drqdbetaf - d(R/q)/d(beta) 
dsfqdbetaf - d(Y/q)/d(beta) 

dalpfdxdotf - d(alpha fitse)/dXdot 
dalpfdzdotf • d(alpha fuse)/dZdot 
dbetadydotf - d(beta)/dYdot 
detamfdxdotf - d(eta)M/dXdot 
detamfdzdotf - d(eta)M/dZdot 
dgamdzdotf - d(gamma)/dZdot 
dmdxdotf - dM/dXdot 
dmdzdotf - dM/dZdot 
dndydotf - dN/dYdot 
drdydotf - dR/dRdot 
dxdxdotf - dX/dXdot 
dxdzdotf - dX/dZdot 
dydydotf - dY/dYdot 
dzdxdotf - dZ/dXdot 
dxdzdotf - dZ/dZdot 

HOVER -- BASIC MR Derivatives 

dalda - dals/dAl 
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daldb - dals/dBl 

daldlamp - dals/d(lambda’) 

daldmu - dals/dmu 

daldp - dals/dp 

daldq - dals/dq 

daldtheto - dals/dtheta(o) 

dblda - dbls/dAl 

dbldb - dbls/dBI 

dbldlamp - dbls/d(lambda’) 

dbldmu - dbls/dmu 

dbldp - dbls/dp 

dbldq - dbis/dq 

dbldtheto - dbls/dtheta(o) 

dbetdydot - d(Bcta)/dYdot 

dchsigda - d(Ch/sigma)/dals 

dchsigdb - d(Ch/sigma)/dbls 

dchsigdlamp - d(Ch/sigma)/d(lambda’) 

dchsigdmu - d(Ch/sigma)/dmu 

dchsigdUieto - d(Ch/sigma)/dtheta(o) 

dctsigdlamp - d(Ct/sigma)/d(lambd; , 

dctsigdmu - d(Ct/sigma)/dmu 

dctsigdtheto - d(Ct/sigma)/dtheta(o) 

dcqsigdlamp - d(Cq/sigma)/d(lambda’) 

dcqsigdmu - d(Cq/sigma)/dmu 

dcqsigdtheto - d(Cq/sigma)/dtheta(o) 

dcysigdb - d(Cy/sigma)/dbls 

dlampdxdot - d(lambda’)/dXdot 

dlampdydot - d(lambda’)/dYdot 

dlampdzdot - d(lambda’)/dZdot 

dmdalsm - dM/dals - main rotor stiffness 

dmudxdot - dmu/dxdot 

drdblsm - dR/dBls - main rotor stiffness 

RIGGING -- Cockpit stick rigging gains 

dalmddela - dals/d(delta aileron) 

dblmddele - dbls/d(delta elev) 

ddelvddelp - rudder deflection per pedal travel 

dphinddelp - dphifnotar thruster sleeve angle)/d(pedal) 

dthetodc * dtheta(0.7)/d(collective) 

dthetomddelc - dthetaOm/d(collective) 

dthetotddelp - dtheta(0.7,tail rotor)/d(pedal travel) 
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APPENDIX E. PLOT LISTING. 

Plots are saved under the following filenames: 

Longitudinal roots - rootlon.met 
Lateral/Directional roots - rootlat.met 
Coupled roots • rootcoup.met 

Longitudinal Cyclic 

Longitudinal Cyclic to U, Hover - cbe2uh.met 
Longitudinal Cyclic to Theta, Hover - cbe2theh.met 
Longitudinal Cyclic to Pitch Rate, Hover - cbe2qh.met 
Longitudinal Cyclic to W, Hover - cbe2wh.met 

Collective 

Collective to U, Hover - ebc2uh.met 
Collective to Pitch, Hover - cbc2theh.met 
Collective to Pitch Rate, Hover - cbc2qh.met 
Collective to W, Hover - cbc2wh.met 

Lateral cyclic 

Lateral Cyclic to Bank, Hover - cba2phih.met 
Lateral Cyclic to Sideslip (v). Hover - cba2vh.met 
Lateral Cyclic to Roll Rate, Hover - cba2ph.met 
Lateral Cyclic to Yaw Rate, Hover - cba2rh.met 
Lateral Cyclic to Yaw, Hover - cba2yh.met 

Pedals 

Pedals to Bank, Hover - cbp2phih.met 
Pedals to Sideslip (v), Hover - cbp2vh.met 
Pedals to Roll Rate, Hover - cbp2ph.met 
Pedals to Yaw Rate, Hover - cbp2rh.met 
Pedals to Yaw, Hover - cbp2yh.met 
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APPENDIX F. SAMPLE OUTPUT. 


Sample output file for the Prouty helicopter at a hover. 

*** RESULTS *** 
prouth 

*** INPUT DATA *** 


Flight Conditions 


Forward velocity = 

0 

kts 

Temperature = 

90 

degs F 

Pressure altitude = 

0 

ft 

Auxiliary thrust = 

0 

lbs 

Fuselage 


Gross weight = 

20000 

lbs 

Equivalent flat plate area = 

19.3 

ft A 2 

Vertical projected area = 

380.0 

ft A 2 

CG height above waterline = 

0.0 

ft 

CG fuselage station = 

0.0 

ft 

CG position rt of buttline = 

0.0 

ft 

IXX = 

35000.C 

) slug ft A 2 

iyy = 

40000.C 

) slug ft A 2 

Izz = 

35000.C 

) slug f t A 2 

Ixz = 

0.0 

slug ft A 2 

Downwash ratio * 

1.50 


Main Rotor 


Number of blades = 

4 


Rotor radius * 

30.0 

ft 

Blade chord = 

2.0 

ft 

Blade twist = 

10.00 

degs 

Blade airfoil = 

HH-02 


Blade lift curve slope = 

5.73 


Blade weight = 

207.0 

lbs 

Rotational velocity = 

21.67 

rads/sec 

Blade grip length = 

4.5 

ft 

Hinge offset = 

1.5 

ft 

Flapping moment of inertia = 

2870.0 

slug ft A 2 

Hub height above waterline = 

7.5 

ft 

Hub fuselage station = 

0.0 

ft 

Hub position rt of buttline = 

0.0 

ft 

Mast incidence & 

0.00 

deg 

Tail rotor (zero if 

NOTAR) 

Number of blades = 

3.0 


Blade chord = 

1.0 

ft 

Blade radius = 

6.5 

ft 

Lift curve slope = 

5.73 


Rotational velocity = 

100.00 

rad/sec 

Flapping moment of inertia = 

- : 

r-lug ft A 2 
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Delta-3 angle « -30.00 deg 

Blade twist > -5.00 deg 

Hub height above waterline « 6.0 ft 

Hub fuselage station » 37.0 ft 

Hub position rt of buttline ■ 0.0 ft 

Wing 

Area - 0.0 ft A 2 

Span * 0.0 ft 

CL . 0.00 

CDo » 0.0000 

Tip cord * 0.0 ft 

Root cord > 0.0 ft 

wing efficiency factor x 0.00 

Zero lift angle x 0.00 deg 

Angle of incidence & 0.00 deg 

Lift curve slope * 0.00 

Height above waterline * 0.0 ft 

Fuselage station x 0.0 ft 

Position right of buttline * 0.0 ft 

Rotor downwash ratio = 0.00 

Fuselage downwash ratio = 0.00 

Horizontal tail 

Area = 18.0 ft A 2 

Span = 9.0 ft 

CL = 0.20 

CDo = 0.0045 

Zero lift angle = 0.00 deg 

Angle of incidence = 0.00 deg 

Lift curve slope * 5.73 

Height above waterline * -1.5 ft 

Fuselage station * 33.0 ft 

Position right of buttline = 0.0 ft 

Dynamic pressure ratio = 0.60 

Rotor downwash ratio = 1.50 

Fuselage downwash ratio = 1.20 

Vertical tail 

Area = 33.0 ft A 2 

Span = 7.7 ft 

CL = 0.20 

CDo = 0.0045 

Height above waterline * 3.0 ft 

Fuselage station & 35.0 ft 

Position right of buttline = 0.0 ft 

Zero lift angle = 0.00 deg 

Maximum Cl = 2.00 

Dynamic pressure ratio = 0.60 

Lift curve slope = 5.70 

Rigging 

Long cyclic pitch/inch defl * 3.33 deg/in 

Lat cyclic pitch/inch defl x 2.05 deg/in 

Collective pitch/inch defl x 1.33 deg/in 

Tail rotor pitch change/defl z -8.70 deg/unit 
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0.09 units 


Max deflection of control 
from neutral for NOTAR * 

Displacement of anti-torque 

control until full rudder * 0.00 units 


*** CALCULATED DATA *** 


Main Rotor 


advance ratio 
inflow parameter wrt TPP 
Tip path angle 
Rotor coning angle 
1st lat cyclic term-Al 
1st long cyclic term-Bl 
lateral flapping 
longitudinal flapping 
Lock number 

Tail Rotor 

tail rotor thrust 
advance ratio 
inf l'- parameter 
Roto; r.ing angle 
lat-, .-1 flapping 
longitudinal flapping 
Lock number 


m 

0.0 


X 

0.000 


X 

0.0 

degs 

s 

7.4 

degs 

X 

0.00 

degs 

= 

-0.00 

degs 

X 

0.0 

degs 

= 

0.0 

degs 

X 

7.2 


(zero if 

NOTAR) 

_ 

1349.5 

lbs 

X 

0.0 


= 

-0.073 


= 

-0.1 

degs 

= 

0.00 

degs 

X 

0.00 

degs 

X 

2.7 



State Matrices 


Longitudinal uncoupled plant (A or F depending on notation) 
States are [u w q theta] 


-0.0106 

0 

0.0038 

0 


0 

-0.2698 

0 

0 


2.^068 

0 

-0.8878 

0.9977 


-32.2000 

0.0002 

0 

0 


Longitudinal uncoupled input matrix (B or G depending on notation) 
Inputs are [longitudinal cyclic, collective, lateral cyclic, pedals] 


1.2248 

-0.0000 

-0.0738 

0 

0 

-5.8773 

0 

0 

0.4338 

0.0000 

0.0261 

-0.0238 

0 

0 

0 

0 


Lateral/directional uncoupled plant (A or F depending on notation) 
States are [v p phi r psi] 


-0.0305 

-0.0064 

0 

0.0130 

0 


-2.6260 

-1.0274 

1.0000 

0.0783 

0 


32.1267 

0 

0 

0 

0 


0.7353 

0.0783 

- 0.0000 

-3.3360 

1.0000 


0 

0 

0 

0 

0 
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Lateral/directional uncoupled input matrix 
on notation) 

Inputs are [longitudinal cyclic, collective 


0.1200 0.0000 

0.0486 0.0000 

0 0 

0 0.3001 

0 0 


0.7533 -1.9736 

0.3049 -0.2101 

0 0 

0 1.2959 

0 0 


(B or G depending 
, lateral cyclic. 


pedals] 


Coupled plant (A or F depending on notation) 
States are [u w q theta v p phi r psi] 


Columns 1 
-0.0106 

through 7 

0 

2.5068 

-32.2000 

-0.0057 

-0.7469 

0 

0 

-0.2698 

0 

0.0002 

0 

0 

2.1710 

0.0038 

0 

-0.8878 

0 

0.0013 

0.2645 

0 

0 

0 

0.9977 

0 

0 

0 

0 

0.0057 

0 

-0.7469 

-0.0000 

-0.0305 

-2.6260 

32.1267 

0.0023 

0 

-0.3023 

0 

-0.0064 

-1.0274 

0 

0 

0 

0.0000 

0 

0 

1.0000 

0 

0 

-0.0198 

0 

0 

0.0130 

0.0783 

0 

0 

0 

0 

0 

0 

0 

0 


Columns 8 through 9 


0 0 

0 0 

0.0266 0 

0.0674 0 

0.7353 0 

0.0783 0 

- 0.0000 0 

-3.3360 0 

1.0000 0 


Coupled input matrix (B or G depending on notation) 

Inputs are [longitudinal cyclic, collective, lateral cyclic, pedals] 


1.2248 -0.0000 

-0.0738 

0 

0 -5.8773 

0 

0 

-0.4338 0.0000 

0.0261 

-0.0238 

0 0 

0 

0 

0.1200 0.0000 

0.7533 

-1.9736 

0.0486 0.0000 

0.3049 

-0.2101 

0 0 

0 

0 

0 0.3001 

0 

1.2959 

0 0 

Uncoupled 

0 

Eigenvalue 

0 

Longitudinal plant 

Root 

wn 

damping 

-1.0155 

1.0155 

1.0000 

0.0585 + 0.3398i 

0.3448 

-0.1698 

0.0585 - 0.3398i 

0.3448 

-0.1698 

-0.2698 

0.2698 

1.0000 
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Lateral/Directional plant 


Root 

wn 

damping 

0 

0 

- 1.0000 

-3.3433 

3.3433 

1.0000 

-1.1784 

1.1784 

1.0000 

0.0639 + 0.40251 

0.4075 

-0.1568 

0.0639 - 0.4025i 

0.4075 

-0.1568 


Coupled Plant 



Root 

wn 

damping 

0 



0 

- 1.0000 

-3.3432 



3.3432 

1.0000 

-1.0970 

♦ 

0.2583i 

1.1270 

0.9734 

-1.0970 

- 

0.2583i 

1.1270 

0.9734 

0.0760 

+ 

0.3983i 

0.4054 

-0.1873 

0.0760 

- 

0.3983i 

0.4054 

-0.1873 

0.0461 

♦ 

0.3582i 

0.3611 

-0.1276 

0.0461 

- 

0.35821 

0.3611 

-0.1276 

-0.2689 



0.2689 

1.0000 







APPENDIX G. SAMPLE OPEN LOOP PLOTS. 


Saa$>le open loop transfer plots 


Open loop response LongitixJmol Cyclic to U. Hov«r 



Open loop response lonqitudinol Cyclic to Theta. Hover 



Angular Frequency (rod/sec) 
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Gain (d8) Gain (OB) 


Open loop response Lot oral Cyclic to Roll Rate. Hover 



Angular Frequency (rod/sec) 


Open loop response Loteral Cyclic to Yaw Rote. Hover 
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Open loop response loterol Cyclic to fa*. Hover 



Angular Frequency (rod/sec) 


Open loop response Pedals to Bonl«, Hover 



Angular Frequency (rod/sec) 
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Coin (<J8) Coin (dB) 



Open loop response Pedals lo Roll Rote, Hover 



10-2 10 -' 10 « 10 ’ 102 
Angular Frequency (rod/sec) 
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Op«n loop response PeOola to Yaw Rata, Hover 



Angular Frequency (rod/sec) 


Open loop rasponsa Padols to Yaw, Hover 



Angular Frequency (rod/sec) 
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Gain (dB) 


Open loop rMponee Longitudinal Cyclic to U. Cruise 



Angular Frequency (rod/sec) 


Open loop response longitudinal Cyclic to Theto. Cruise 



Angular Frequency (rod/sec) 
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Open loop response Longitudinol Cyclic to Pitch Rote. Cruise 



Angular Frequency (rod/see) 


Open loop response Longitudinal Cyclic to W. Cruise 



Angular Frequency (rod/sec) 
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Gain (dB) Gain (dB) 


Open loop response Loterol Cyclic to Bonk, Cruise 



Open loop response laterol Cyclic to Sideslip (v). Cruise 
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Cain 


Open loop response lateral Cyclic to Ro* Roto. Cruise 



Angular Frequency (rod/sec) 


Open loop response Lateral Cyclic to Yqw Rote, Cruise 



Angular Frequency (rod/sec) 
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Open loop response laterol Cyclic to You. CruiM 



Angular Frequency (rod/sec) 


Open loop response Pedate to Bonk, Cruise 



Angular Frequency (rod/sec) 
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APPENDIX H. COMPUTER CODE. 


Conqputer code for JANRAD Stability 

A. STAB.N 
% STAB.M 

% NPS Helo Preliminary Design Program 
% Stability and Control Routines 
% Written by MAJ Walter M. Wirth, Jr. 

% September 1993 

% This program was designed as an interactive preliminary 
% design tool for stability and control analysis of a single 
% main rotor conventional or compound helicopter. The 
% program provides stability derivatives, roots of the aircraft 
% plant, plots of various control parameters as well as open loop 

% control bandwidths from control inputs to aircraft reponse. 

clear 
load temp 

eval(['load ', filenamel]); 
clc 

disp(' ') 
disp(' ') 

disp(' *** STABILITY AND CONTROL ROUTINE ***') 

disp(' ') 
disp(' ') 
pause(3) 

%%%%%%%*%%%%%%%%%%%%%%%%%%*%%%*%%%%%%%%%%%%%%»%»%»%%%%%%%%%%%%%%% 

% *** If editing an existing file: get file name, display edit 

% menu, allow changes to selected variables, and save under 

% desired file name. Loads to and saves from current 

% directory as a .mat file. *** 

if exist('lb '); 

answerOsl; 
else 

answer0=2; 
end 

if answer0==l, 
check2=l; 
while check2 > 0 
clc 
disp(' ') 

disp(' *** STABILITY AND CONTROL MENU ***') 

disp(' *** ADDITIONAL PARAMETERS (1 of 3) ***') 

disp(' ') 

disp(' Main Rotor') 

disp(' 1. flapping mean of inertia 2. hub height above waterline') 

disp(' 3. hub fuselage station 4. hub posn right of buttline') 

disp(' 5. mast incidence') 

% 

disp (' ') 

disp(' Tail Rotor (enter zeros (0) if using NOTAR)') 
disp(' 6. height above waterline 7. hub fuselage station') 

disp(' 8. posn right of buttline 9. number of blades') 
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disp ('10. blade chord 
disp('12. lift curve slope 
disp('l*. flap morn of inertia 
disp('16. blade twist') 

% 


ll. blade radius’) 

13. rotational velocity') 
15. delta-3 angle') 


diep(' ') 

disp(' Verticle Pin') 

disp('17. height above waterline 18. fuselage station') 
disp('19. posn right of buttline 20. alpha zero lift’) 
disp('21. CL> max 22. dynamic pressure ratio') 

disp('23. lift curve slope') 

% 

disp(' ') 

disp('0. NO CHANGES') 

choice-input('Input the parameter to change: '); 
it choicec«l, 
clc 

disp(' ') 
lb 


tempi-lb; 

Ib-input('Blade flapping moment of inertia (slug ft A 2): ') 
if isempty(Ib), 

Ib=templ; 

end 

clear tempi 
elseif choice==2, 
clc 

disp(' ') 
hmd 

tempishmd; 

hmd=input('Hub height above reference datum/waterline (ft) 
if isempty(hmd), 
hmd=templ; 

end 

clear tempi 
elseif choice==3, 
clc 

dispC ') 
lmd 

tempi-lmd; 

lmd=input('Hub fuselage station (ft): '); 
if isempty(lmd), 
lmd-templ; 

end 

clear tempi 
elseif choice=»4, 
clc 

disp(' ') 
ymd 

tempi-ymd; 

ymd=input('Hub position right of buttline (ft): '); 
if isempty(ymd), 
ymd-templ; 
end 

clear tempi 
elseif choice-*5, 
clc 

disp(' ') 
im*57.3 

templ*im*57.3; 
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im*input('Mast incidence (negative forward - deg): ')/57.3; 
if isempty(im), 
im*templ/57.3; 

end 

clear tempi 
elseif choice-*6, 
clc 

disp(' ') 
htd 

tempi*htd; 

htd*input ('Tail rotor height above reference datum/waterline (ft) 

if isenpty (htd) , 
htd*templ; 
end 

clear tempi 
elseif choice**7, 
clc 

disp(' ') 
ltd 

tempi*ltd; 

ltd*input('Tail rotor fuselage station (ft): '); 
if isenpty(ltd), 
ltdstempl; 

end 

clear tenpl 
elseif choice=*8, 
clc 

disp(' ') 
ytd 

tenpl=ytd; 

ytd=input('Tail rotor position right of buttline (ft): ’); 
if isenpty(ytd) , 
ytd=tempi; 
end 

clear tempi 
elseif choice=«9, 
clc 

disp(' ') 
bt 

templ*bt; 

bt=input('Number of tail rotor blades: '); 
if isempty(bt), 
bt*templ; 

end 

clear tenpl 
elseif choiceaalO, 
clc 

disp(' ') 
cot 

tempi*cot; 

cot*input('Blade chord (ft): '); 
if isenpty(cot), 
cot*templ; 

end 

clear tempi 
elseif choice**ll, 
clc 

disp(' ') 

Rt 
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tempi»Rt; 

Rt-input('Tail rotor blade radius (ft): '); 
if isempty(Rt), 

Rtatempl; 

end 

clear tempi 
elseif choice««l2, 
clc 

disp(* ') 

at 

tempisat; 

at«input('Average lift curve slope of tail rotor: '); 
if isempty(at), 
atatempl; 

end 

clear tempi 
elseif choices*13, 
clc 

disp(' ') 
ohmt 

tempi=ohmt; 

ohmt■input('Rotational velocity of tail rotor (rad/sec): '); 
if isempty(ohmt), 
ohmt=tempi; 

end 

clear tempi 
elseif choice==14, 
clc 

disp(' ') 

Ibt 

templ=Ibt; 

Ibt=input('Blade flapping moment of inertia (slug ft A 2) : ') ; 
if isempty(Ibt), 

Ibt=templ; 

end 

clear tempi 
elseif choice**15, 
clc 

disp(' ') 
delta3*57.3 
templ®delta3*57.3; 

delta3»input('Delta-3 angle (deg): 'J/57.3; 
if isempty(delta3), 
delta3=templ/57.3; 

end 

clear tempi 
elseif choice»16, 
clc 

disp(' ') 
thetalt*57.3 
templathetalt*57.3; 

thetaltsinput('Blade twist (deg): 'J/57.3; 
if isempty(thetalt), 
thetaltatempl/57.3; 

end 

clear tempi 
elseif choice«l7, 
clc 

disp(' ') 
hvd 
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tempi*hvd; 

hvd-input('Height above reference datum/waterline (ft): '); 
if isempty(hvd), 
hvd-tempi; 

end 

clear tempi 
elaeif choice»«18, 
clc 

disp(* *) 
lvd 

tempi-lvd; 

lvd-input('Fuselage station (ft): '); 
if isempty(lvd), 
lvd-templ; 

end 

clear tentpl 
elseif choice-«19, 
clc 

disp(' ') 
yvd 

templ-yvd; 

yvd-input('Position right of buttline (ft): '); 
if i6empty(yvd) , 
yvd-f ""npl ; 

end 

clear tempi 
elseif choice-=20, 
clc 

dispC ') 
alplov*57.3 
templ-alplov*57.3; 

alplov-input('Zero lift angle for vertical tail (deg): ')/57.3 
if isempty(alplov), 
alplov-templ/57.3; 

end 

clear tempi 
elseif choice--21, 
clc 

disp(' ') 

clvertmax 

tempi-clvertmax; 

clvertmax-input('Maximum Cl for vertical tail: '); 
if isempty(clvertmax) , 
clvertmax-templ 

end 

clear tempi 
elseif choice«*22, 
clc 

disp(' ') 

qvq 

tempi-qvq; 

qvq-input('Dynamic pressure ratio (pg 489 Prouty): '); 
if isempty(qvq) , 
qvq-tempi; 

end 

clear tempi 
elseif choice--23, 
clc 

disp(' ') 

av 
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tampl«av; 

av* input ('Lift curve elope of vertical tail: '); 
if isempty(av), 
av.templ; 

end 

clear tenipl 
elseif choice««0, 
check2»0 
clc 

else 

disp(' ') 

diep('enter a displayed number ...press any key to continue') 

pause 

end 

end 

t 


check2«l; 


while 
clc 
disp(' ') 
disp(' 

check2 > 0 

*** STABILITY AND CONTROL MENU ***') 

disp(' 

*** ADDITIONAL PARAMETERS (2 of 3) ***') 

disp(' ') 

disp(' Horizontal Tail') 
disp(* 1. height above waterline 

2. 

fuselage station') 

disp(' 3. 

posn right of buttline 

4. 

alpha zero lift') 

disp(' 5. 

angle of incidence 

6. 

lift curve slope') 

disp(' 7. 

dynamic pressure ratio 

8. 

rotor downwash ratio') 

disp(' 9. 

fuselage downwash ratio' 

) 


disp(' ') 
dispC wing') 

disp!'10. height above waterline 

11. 

fuselage station') 

disp('12. 

posn right of buttline 

13. 

alpha zero lift') 

disp ('14 . 

angle of incidence 

15. 

lift curve slope') 

disp('16. 

tip cord 

17. 

root cord') 

disp('18. 

ft. 

rotor downwash ratio 

19. 

fuselage downwash ratio 

disp(' ') 

disp('0. NO CHANGES') 




choice=input('Input the parameter to change: '); 


if choice«l 
clc 

disp(' ') 
hhd 

tenq?l«hhd; 

hhd=input('Height above reference datum/waterline (ft): '); 
if isempty (hhd) , 
hhd-templ; 
end 

clear tempi 
elseif choice«*2, 
clc 

disp(' ') 
lhd 

tempi■lhd; 

lhd«input('Fuselage station (ft): '); 
if isempty(lhd), 
lhdatempl; 
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end 

clear tempi 
elseif choice..3, 
clc 

disp(' ') 
yhd 

tempi.yhd; 

yhd*input('Position right of buttlins: ') ; 
if isempty(yhd), 
yhd.templ; 
end 

clear tempi 
elseif choice**4, 
clc 

diep(' ') 
alploh*57.3 
tempi.alploh*57.3; 

alploh.input('Zero lift angle for horizontal tail (deg): ')/57.3; 
if isempty(alploh), 
alploh.templ/57.3; 

end 

clear tempi 
elseif choice.*5, 
clc 

disp(' ') 
ih*57.3 

templ.ih*57.3; 

ih*input('Angle of incidence of horizontal tail (deg): ')/57.3; 
if isempty(ih), 
ih»templ/57.3; 

end 

clear tempi 
elseif choice==6, 
clc 

disp(' ') 
ah 

tempi*ah; 

ah*input('Lift curve slope of horizontal tail: '); 
if isempty(ah), 
ah.templ; 

end 

clear tempi 
elseif choice=*7, 
clc 

disp(' ') 
qhq 

tempi*qhq; 

qhq*input('Dynamic pressure ratio (pg 489 Prouty): '); 
if isempty(qhq), 
qhq*templ; 

end 

clear tempi 
elseif choice*«8, 
clc 

disp(' ') 
vhvl 

tempi*vhvi; 

vhvl.input('Rotor downwash ratio (pg 489 Prouty): '); 
if isempty(vhvl), 
vhvl.tempi; 
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•nd 

clear tempi 
elaeif choice»«9, 
clc 

disp <' ') 

detafdalpfh 

tempi-detafdalpfh; 

detafdalpfh-input('Fuselage downwash ratio (pg 489 Prouty): '); 
if isempty(detafdalpfh), 
detafdalpfh-tempi; 

end 

clear tenqpl 
elseif choice«*10 
clc 

disp(' ') 
hwd 

templ-hwd; 

hwd-input('Height above reference datum/waterline (ft): '); 
if isenpty(hwd) , 
hwd=tempi; 

end 

clear tempi 
elseif choice==ll 
clc 

disp(' ') 
lwd 

tempi-lwd; 

lwd=input('Fuselage station (ft): '); 
if isempty(lwd), 
lwd=templ; 

end 

clear tempi 
elseif choice-*12 
clc 

dispC ') 
ywd 

tempi-ywd; 

ywdsinput('Position right of buttline (ft): '); 
if isempty(ywd), 
ywd=templ; 
end 

clear tempi 
elseif choice=*13, 
clc 

disp(' ') 
alplow*57.3 
templ*alplow*57 .? • 

alplow-input('Zer^ xift angle for wing (deg): *)/57.3; 
if isempty(alplow), 
alplow-templ/57.3; 

end 

clear tempi 
elseif choice--l4, 
clc 

disp{* ') 
iw*57.3 

templ-iw*57.3; 

iw-input('Angle of incidence of wing (deg): ')/57.3; 
if isempty(iw), 
iw«templ/57.3; 
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•nd 

clear tempi 
elseif choice**15, 
clc 

diep(' ') 

aw 

tempi*aw; 

aw-input('Lift curve elope of wing: '); 
if isempty(aw), 
aw*tempi; 

end 

clear tespl 
elaeif choice-*l6, 
clc 

diep(' ') 
ctw 

tempi*ctw; 

ctw*input('Tip cord (ft): '); 
if isenqpty (ctw), 
ctw*tempi; 

end 

clear tempi 
elseif choice**17, 
clc 

disp(' ') 
erw 

tempi*crw; 

erw*input('Root cord (ft): '); 
if isempty(erw), 
crw*templ; 

end 

clear tempi 
elseif choice**18, 
clc 

disp(' ') 
vwvl 

templevwvl ; 

vwvi*input('Rotor downwash ratio (pg 489 Prouty): '); 
if isempty(vwvl), 
vwvl*tempi; 

end 

clear tempi 
elseif choice**19, 
clc 

disp(' ') 

detafdalpfw 

templ*detafdalpfw; 

detafdalpfw*input('Fuselage downwash ratio (pg 489 Prouty): ') 
if isempty(detafdalpfw), 
detaf dalpf watempl 

end 

clear tempi 
elseif choice*«0, 
check2*0 
clc 

else 

dispC ') 

disp('enter a displayed number ... press any key to continue') 

pause 

end 







•nd 

% 

check2*l; 
while check2 > 0 
clc 
disp(* ') 

disp(' *** STABILITY AND CONTROL MENU ***') 

dispC *** ADDITIONAL PARAMETERS (3 of 3) ***') 

% 

disp(' ') 

dispC CG location and Inertias/fuselage parameters') 
dispC 1. eg ht. above waterline 2. eg fuselage station') 
dispC 3. eg posn rt of buttline 4. Ixx') 
dispC 5. Iyy 6. Izz') 

dispC 7. Ixz 8. fuselage downwash ratio') 

» 

disp(' ') 

dispC NOTAR if available (enter zeros if using tail rotor)') 
disp(' 9. height above waterline 10. boom fuselage station') 
disp('ll. boom position left ref 12. NOTAR diameter ') 
disp('13. swirl angle at boom 14. NOTAR max force') 
disp('IS. thruster fuselage station') 

% 

disp(' ') 
disp(' Rigging') 

disp('16. B1 main/in defl (del e) 17. A1 main/in defl (dela)') 

disp('18. thetaOm/in defl (del c) 19. thetaOt/pedal defl (del r or p)') 

disp('20. NOTAR sleeve twist/defl 21. max rudder defl') 

dispC ') 

dispC 0. NO CHANGES') 

choice*inputCInput the parameter to change: '); 
if choice**l, 
clc 

dispC ') 
zcg 

tempi*zcg; 

zcg=input('CG height above reference datum/waterline (ft): '); 
if isempty(zcg), 
zcg-templ; 

end 

clear tempi 
elseif choice**2, 
clc 

dispC ') 
xcg 

tempi*xcg; 

xcg*input('CG Fuselage station (ft): '); 
if isempty(xcg), 
xcg-templ; 
end 

clear tempi 
elseif choice**3, 
clc 

dispC ') 
yeg 

tempi*ycg; 

yeg*input('CG position right of buttline (ft): '); 
if isempty(yeg), 
ycgatempl; 
end 
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clear tempi 
eleeif choice**4, 
clc 

dispC ') 

Zxx 

templ»Ixx; 

Ixx*inputCIxx (slug ft*2): '); 
if isempty(Ixx), 

Ixx*tempi; 

end 

clear tempi 
eleeif choice*«5, 
clc 

disp(' ') 
lyy 

templ*Iyy; 

lyy■input('lyy (slug ft A 2): '); 
if isempty(lyy), 

Iyy.templ; 

end 

clear tempi 
elseif choice**6, 
clc 

dispC ') 

IZ2 

templ*Izz; 

Izz*input('Izz (slug ft“2) : 
if isempty(Izz), 

Izz*templ; 

end 

clear tempi 
elseif choice*=7, 
clc 

disp(' ') 

Ixz 

tempi* 1 x 2 ; 

Ixz*inputCIxz (slug ft*2): '); 
if isempty(Ixz), 

Ixz*tempi; 

end 

clear tempi 
elseif choice**8, 
clc 

dispC ') 
vfvl 

tempi*vfvl; 

vfvi*input CDownwash ratio for fuselage (page 513 Prouty): '); 
if isempty(vfvl), 
vfvl*templ; 
end 

clear tempi 
elseif choice**9, 
clc 

disp(' ') 
htnd 

tempi*htnd; 

htnd*inputC Height above reference datum/waterline (ft): '); 
if isempty(htnd), 
htnd*templ; 

end 
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clear tempi 
elseif choice**lO, 
clc 

disp(' ') 
lend 

tempi*ltnd; 

ltnd*input('Fuselage station (ft): '); 
if isempty(ltnd), 
ltnd-tempi; 

•nd 

clear tempi 
elseif choice**ll, 
clc 

disp(' ') 
ytnd 

tempi-ytnd; 

ytnd*input('Position right of buttline (ft): '); 
if isempty(ytnd), 
ytnd-templ; 

end 

clear tempi 
elseif choice--12, 
clc 

disp(' ') 
dian 

tempi-dian; 

dian*input('NOTAR boom diameter (ft): '); 
if isempty(dian), 
dian*templ; 

end 

clear tempi 
elseif choice--13, 
clc 

disp (' ') 
swirl*57.3 
templ-swirl*57.3; 

swirl-input('Swirl angle at boom(deg): ')/57.3; 
if isempty(swirl), 
swirl*templ/57.3; 

end 

clear tempi 
elseif choice-*14, 
clc 

disp(' ') 

Ytmaxn 

tempi-Ytmaxn; 

Ytmaxn*input('Maximum NOTAR thruster force (lbs): '); 
if isempty(Ytmaxn), 

Ytmaxn*tempi; 

end 

clear tempi 
elseif choice**l5, 
clc 

disp(' ') 
lttnd 

templ-lttnd; 

lttnd*input('Thruster fuselage station (ft): '); 
if isempty(lttnd), 
lttnd-templ; 

end 
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clear tenpl 
•laaif choice-*16, 
clc 

disp(' ') 
dblmddele*57.3 
templ-dblmddele*57.3; 

dblmddele*input('Long cyclic pitch per inch defl (deg/in): 

')/57.3; 

if ieempty(dblmddele) , 
dblmddele«templ/57.3; 

end 

clear tempi 
elseif choice««l7, 
clc 

diep(' ') 
dalmddela*S7.3 
tempi-dalmddela*57.3; 

dalmddela*input('Lateral cyclic pitch per inch defl (deg/in): 

')/57.3; 

if ieempty(dalmddela), 
dalmddela*templ/57.3 ; 

end 

clear tempi 
elseif choice»l8, 
clc 

disp(' ') 
dthetomddelc*57.3 
templ»dthetomddelc»57.3; 

dthetomddelc*input('Collective pitch per inch defl (deg/in): 

')/57.3; 

if isempty(dthetomddelc), 
dthetomddelc*templ/57.3; 

end 

clear tempi 
elseif cnoicesKl9, 
clc 

disp <' ') 
dthetotddelp*57.3 
templ*dthetotddelp*57.3; 

disp('Tail rotor pitch change per inch defl or percentage of 

twist') 

disp('Enter 0 (zero) if using NOTAR') 

dthetotddelpxinput('(deg/in or deg/deg of twist): 'J/57.3; 
if isempty(dthetotddelp), 
dthetotddelpxtempl/57.3; 

end 

clear tempi 
elseif choice**20, 
clc 

disp(' ') 
sidearm/2 
templxsidearm/2; 

disp('Maximum deflection of <mti-torque from neutral for NOTAR, 

enter') 

sidearm»input('1000 if using tail rotor (deg or inch travel): 

') *2 ; 

if isempty(Bidearm), 
sidearmxtempl*2; 
if sidearm>«0,sidearm«le3,end 
end 
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clear tempi 

dphinddelp-pi/sidearm; %% pi rad sleeve twist/sidsarm defl 
elseif choice«»2l, 
clc 

disp(' ') 
maxr 

tempi-maxr; 

disp{'Displacement of anti-torque control until full rudder') 
disp(' deflection. Enter 0 (zero) if rudder is fixed') 
maxr-input(' (deg or inch travel): '); 
if isempty(maxr), 
maxr-tempi; 
end 

clear tempi 
elseif choice=»0, 
clc 

check2«0; 
else 

disp(' ') 

disp('enter a displayed number ... press any key to continue') 
pause 
end 
end 

% 

dispC ') 
disp (' ') 

disp(' *** SAVE INSTRUCTIONS **♦*) 

disp(' ') 

disp('A. Save the new data to a specified file name.') 

disp('B. Do not use an extension or quotations.') 

disp('C. Use letter/number combinations of 6 characters or less.') 

disp('D. The file will be saved with a ".mat" extension.’) 

disp(' ') 

dispf'E. If you made no changes or want the same name, press enter.') 

disp (' ') 

disp('ex: desig2') 

filename-filenamel; 

filenamel-input('save file as: ','s'); 
if isempty(filenamel), 
filenamelsfilename; 

end 

clear check 

evalU'save '.filenamel]); 
check*0; 

% end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*%«%%%%%%%%%%%%%%%%%%%%%%% 

% *** If creating a new file: get input for required variables 

% and save under desired file name. Saves to current 

% directory as a .mat file. *** 


else 

check4-l; 
while check4>0; 
clc 

disp('Do you want to use a Tail Rotor or NOTAR?') 
tenp«input('Tail Rotor«0, NOTAR-1: '); 
if temp«-0 
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notar-O; 
checks-0; 

•laaif temp«»l 
notar-1; 
checks-0; 

•laa 

disp(' ') 

disp ('Enter a 0 or 1') 

disp('press any key to continue...') 

pause 

end 

end 

checks-l; 
while check4>0; 
clc 

disp('Do you want to use a controlable vertical tail?') 
temp-input (' No-0, Yes-1: ') ; 
if temp—0 
ctail-0; 
check.4-0 ; 
elseif temp--l 
ctail-1; 
check4*0; 
else 

disp(' ') 

disp('Enter a 0 or 1') 
disp('press any key to continue...') 
pause 
end 
end 

if Swing<.1 
wing-0; 
else 

wing-1; 
end 
clc 

disp('Main rotor') 
disp (' ') 

lb-input('Blade flapping moment of inertia (slug ft*2) : ’); 
hmd-input('Hub height above reference datum/waterline (ft): '); 
lmd-input('Hub fuselage station (ft): '); 
ymd*input('Hub position right of buttline (ft): '); 
im-input('Mast incidence (negative forward - deg): 'J/57.3; 

clc 

if notar--0 

disp('Tail rotor') 
disp(' ') 

htd-input(’Tail rotor height above reference datum/waterline (ft) : 

'); 

ltd-input('Tail rotor fuselage station (ft): '); 
ytd-input('Tail rotor position right of buttline (ft): '); 
bt-input('Number of tail rotor blades: '); 
cot-input('Tail rotor blade chord (ft): '); 

Rt«input('Tail rotor blade radius (ft): '); 
at-input('Average lift curve slope of tail rotor: '); 
ohmt-input('Rotational velocity of tail rotor (rad/sec): '); 
Ibt-input('Tail rotor blade flapping moment of inertia (slug 
ft*2): '); 

delta3-input('Delta-3 angle (deg): ')/57.3; 
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thetalt*input('Blade twist (deg): ')/57.3; 
htnd*0;ltnd*0;ytnd*0;dian*0;swirl>0; Ytmaxn*0;lttnd*0; 
elseif notar*«l 
clc 

disp('NOTAR') 
disp(' ') 

htnd*input('Height above reference datum/waterline (ft): '); 

ltnd>input('Fuselage station (ft): '); 

ytnd*input('Position right of buttline (ft): '); 

dian*input('NOTAR boom diameter (ft): '); 

swirl*input('Swirl angle at boom (deg): ')/57.3; 

Ytmaxn*input('Maximum NOTAR thruster force (lbs): '); 
lttndsinput('Thruster fuselage stationerence (ft): '); 
htd* 0;1td* 0;ytd=0;bt*0;cot * 0;Rt * 0;at * 0;ohmt*0; 
lbt*0;delta3«0;thetalt*0; 

end 

clc 

disp('Verticle tail') 
disp(' ') 

hvd*input('Height above reference datum/waterline (ft): '); 

lvd*input('Fuselage station (ft): '); 

yvdainput('Position right of buttline (ft): '); 

alplovfiinput('Zero lift angle for vertical tail (deg): ')/57.3; 

clvertmax=input('Maximum Cl for vertical tail: '); 

qvq=input('Dynamic pressure ratio for tail (pg 489 Prouty): '); 

av*input('Lift curve slope of vertical tail: '); 

clc 

disp('Horizontal tail') 
disp(' ') 

hhd«input('Height above reference datum/waterline (ft): '); 
lhd=input('Fuselage station (ft): '); 
yhdainput('Position right of buttline: '); 

alploh^input('Zero lift angle for horizontal tail (deg): ')/57.3; 
ih*input('Angle of incidence cf horizontal tail (deg): ')/57.3; 
ah*input('Lift curve slope of horizontal tail: '); 
qhq=input('Dynamic pressure ratio for tail (pg 489 Prouty): '); 
vhvl*input('Rotor downwash ratio for h-tail(pg 489 Prouty): '); 
detafdalpfh*input('Fuselage downwash ratio for h-tail (pg 489 
Prouty): '); 
if wing*=l 
clc 

disp(’Wing') 
dispC ') 

hwd=input('Height above reference datum/waterline (ft): '); 

lwd*input('Fuselage station (ft): '); 

ywd*input('Position right of buttline (ft): '); 

alplow*input('Zero lift angle for wing (deg): ')/57.3; 

iw*input('Angle of incidence of wing (deg): ')/57.3; 

aw*input('Lift curve slope of wing: '); 

ctw*input('Tip cord (ft): '); 

erw*input('Root cord (ft): '); 

vwvl*input('Rotor downwash ratio for wing (pg 489 Prouty): '); 
detafdalpfw*input('Fuselage downwash ratio for wing (pg 489 
Prouty): '); 
elseif wing**0 

hwd*0 ; lwd*0;ywd*0;alplow*0;iw*0;aw*0 ; Ctw-0;crw*0; 
vwvl*0;detafdalpfw*0; 

end 

clc 

disp('CG location') 
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disp(' ') 

zcg-input ('CG height above reference datuin/waterline (ft): '); 

xcg-input('CG Fuselage station (ft): '); 

ycg«input('CG position right of buttline (ft): '); 

clc 

disp('Fuselage moments of inertia/downwash parameter') 
disp(' ') 

Ixx-input('Ixx (slug ft A 2): '); 

lyy-input('Iyy (slug ft A 2): ') 

Izzsinput('Is* (slug ft A 2): '); 

Xxz-input{'Ixz (slug ft A 2): '); 

vfvl-input('Downwash ratio for fuselage (page S13 Prouty): '); 

clc 

disp(' Rigging') 
disp(' ') 

dblmddele-input{'Long cyclic pitch per inch defl (deg/in): 

'1/57.3; 

dalmddela-input('Lateral cyclic pitch per inch defl (deg/in): 

')/57.3; 

dthetomddelc>input('Collective pitch per inch defl (deg/in): 

')/S7.3 ; 

if notar**0 
dispC ') 

disp('Tail rotor pitch change per inch defl or percentage of 

twist') 

dthetotddelp=input(' (deg/in or deg/deg of twist): 'J/57.3; 
dphinddelpsO;sidearm=1000; 
elseif notar«*l 
disp(' ') 

disp('Max deflection of anti-torque from neutral for NOTAR, enter 

') 

sidearmsinput(' 1000 if using tail rotor (deg or inch travel): 

') * 2 ; 

if sidearm«=0,sidearm*1000,end 
dphinddelp=pi/sidearm; 

end 

if ctail»i 
disp(' ') 

disp('Displacement of anti-torque control until full rudder') 
maxr*input(' deflection (deg or inch travel): '); 
elseif ctail**0 
maxrsO 

end 

clc 

eval(['save ',filenamel]); 
save stabtemp filenamel 
check-0; 
end 

if diajn--0; 
notar-0; 
if bt»0; 

dispC You must have a tail rotor or NOTAR/thruster! ') 
check-1 
end 
else 

notar-1; 
end 
clc 

disp(' ') 
disp(' ') 
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disp{' *** DATA ENTRY COMPLETE ***') 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

disp<* ') 

disp(' *** EVALUATING STABILITY DERIVATIVES •**') 

pause(3) 

if Vinf<20 

hover % call hover routine 
elseif Vin£>*20 

cruise % call cruise routine 

end 

stabout 

%%»»»%%»%%%%%%%%%%%%%%%%%%%%%%%»%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 





I. TRXM.K 

% TRIM.M 

% Rotor trim subroutine 
disp(' ') 

disp(' *** EXECUTING ROTOR TRIM ROUTINE ***') 

% *** calculation of required parameters *** 

rho».002377*(-.000031*PA+ (-.002*temp+1.118)) ; 

% *** first guess at rotor profile drag ( H force) *** 

if Vinf < 16.9, 

Drotor«0; 

else 

Drotor«Vinf*(rho/.002377); 
end 

q»0.5*rho*Vinf a 2; 

Adisk«pi*R A 2; 

Vtip«omega*R; 

Dfuse»q*Afh; 

CDwing*CDowing+(CLwing A 2/(ewing*pi*(bwing A 2/Swing))); 
CDhoriz*CDohoriz+(CLhoriz A 2/(.8*pi*(bhoriz A 2/Shoriz))); 
CDvert»CDovert+(CLvert A 2/(.8*pi*(bvert A 2/Svert))); 
Dwing*q*CDwing*Swing; 

Dhoriz*q*CDhoriz*Shoriz; 

Dvert*q*CDvert*Svert; 

Oftotals(Dfuse+Dwing+Dhoriz+Dvert)-Taux; 

Lwing*q*CLwing*Swing; 

Lhoriz=q*CLhoriz*Shoriz; 

Lvert«q*CLvert*Svert; 

Lf total-Lwing+Lhoriz; 

alphaT*atan2((Dftotal+Drotor),(GW-Lftotal)); 
mu*Vinf*cos(alphaT)/Vtip; 


% *** thrust calculation *** 

if Vinf < 16.9, 

T=(1+(0.3*Afv/Adisk))*GW; 
CTsT/(Adisk*rho*Vtip A 2); 
else 

T«(GW-Lftotal)/cos(alphaT); 
CT«T/(Adisk*rho*Vtip A 2); 
end 


% *** setup blade radius elements, azimuth elements, 

% induced velocity distributions, and determination 

% of coning angle and tip loss parameter *** 

B«l-(sqrt(2*CT)/b); 

Reff*B*R; 

RbarsReff*e; 

dr« (Ref f-grip) /nbe,- 

r«grip:dr:Reff-dr;,r«r4dr/2; 

rTl«0.7;,% *** first guess at rT *** 

RbarT«rTl*Rbar; 
mblade-wblade/32.17; 
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betao*asin((T/b*RbarT- (.5*(R-e)♦•) *wblade) id.5* (R-e)♦•) A 2*coega A 2*mblad 

•)); 

betat*twist*(0.7-(r/R)); 

psi*0:360/naz:360-360/naz;,psi*psi'/57.3; 


% *** induced velocity determination *** 

if Vinf < 16.9, 

A-4*pi; 

Bv«(b/2)*omega*a*cblade; 

Tv*0; 

delT*T-Tv; 

while abs(delT) > .Ol*T 

thetavmtwist*(0.7-(r/R))+thetao; 

C*(-b/2)*cblade*omega A 2*r*a.*thetav; 
vi*(-Bv+sqrt(Bv A 2-(4*A*C)))/(2*A); 

dTv* (b/2) *rho* ((omega*r) . A 2)*a.* (thetav- (vi. / (omega*r))) *cblade*dr; 
Tv*sum (dTv); 
delT*T-Tv; 
if delT < 0, 

thetao*thetao-0.5*thetao*abs(delT/T); 
else 

thetao*thetao+0.5*thetao*abs (delT/T) ; 
end 

end 

else 

lamdaT*[l-2*mu*sin(alphaT) mu A 2*(sin(alphaT) A 2+l)-2*mu A 3*sin(alphaT) 
mu A 4*sin(alphaT) A 2-CT A 2/4]; 

lamdaT*max(real(roots(lamdaT))); 
vi*lamdaT*Vtip-Vinf*sin(alphaT); 
vi*vi*ones(r); 
end 


% *** first guess at theta *** 

thetalc*0.035*((0.0006e-3*Vinf A 2+0.244e-3*Vinf)/0.105); 
thetals*-0.087*((0.0006e-3*Vinf A 2+0.244e-3*Vinf)/0.105); 
theta*thetao+thetalc.*cos(psi)+thetals.*sin(psi); 


% *** rotor trimming routine *** 

disp(' ') 
disp(' ') 

disp(' *** TRIMMING COLLECTIVE ***') 

k*l; 

errorO*(T*.02)+1; 

while abs(errorO) > T*.02 
Tpsi-zeros(psi); 
thrcalc 

error0*T-(mean(Tpsi)*b); 
if errorO < -T*.02, 

thetao*thetao-0.35*thetao*abs(1.5*errorO/T)*(l-mu); 
elseif errorO > T*.02, 

thetao*thetao+0.35*thetao*abs(1.5*errorO/T)*(l-mu); 
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•nd 

theta*thetao+thetalc.+cos(psi)+thetals.*sin(psi); 
if k > 1, 

if abs(errorO) > abs(errorl), 
clc 

disp(' ') 
disp(' 'i 

disp('This configuration will not trim') 
disp('Try a lower airspeed or a new design') 
disp(' ') 

disp('Program execution terminated*) 
disp(' ') 

error('*** END OF PROGRAM **•') 
end 
end 

errorl*errorO; 
k-k+1; 
end 


disp(' ') 

dispC *** TRIMMING CYCLIC ***') 

tO«clock; 
k*l; 

errorO*(((T/b)*rTl*(R-grip))*.04)+1; 

while errorO > ((T/b)*rTl*(R-grip))*.04 

time-etime(clock,tO); 
if time > IS, 
di8p(' ') 

disp('still trimming ...') 
t0*clock; 
end 

Mpsi(:,k)<zeros(psi); 
tmcalc 

theta*[theta theta(:,k)]; 

Mpsi* [Mpsi Mpsi (:, k)) ; 

% *** calculation of initial dthetadM *** 

if k < 2, 

theta(:,k+1)*theta(:,k)+0.25/57.3; 
l^psi (:, k+1) -zeros (psi) ; 
k*k+l; 
tmcalc 

k.fc .\• 

dthetadM* (theta(;,k+1)-theta(:,k))./(Mpsi(:,k+1)-Mpsi(:,k)) 
end 

% *** calculation of M first harmonic parameters *** 

Mlc*2*sum (*$>si (:, k) . *cos (psi)} /naz ; 

Mls*2*sum(Mpsi(:,k).*sin(psi))/naz; 

% *** removal of first harmonic terms from t^psi *** 

Mpsi(:,k+1)-Mpsi(:,k)-Mlc.*cos(psi)-Mis.*sin(psi); 


120 









delM-Mpsi (:,k+l) -Mpsi (: ,k) ; 
errorO-max(delM)-min(delM); 
if k > 1, 

if errorO > •rrorl, 
clc 

disp(' ') 
disp(' ') 

disp('This configuration will not trim') 
dispi'Try a lower airspeed or a new design') 
disp(' ') 

disp('Program execution terminated') 
disp(' ') 

error{'*** END OF PROGRAM »»*') 
end 
end 

error1-errorO; 

% *** calculation of new theta *** 

delM=0.5*(1-mu)*delM; 

theta(:,k+1)-theta(:,k)+(dthetadM.*delM); 

if errorO <» ((T/b)*rTl*(R-grip))V04, 
thetalc-2*sum(theta(: ,k).*cos(psi))/naz; 
thetals-2*sum(theta(:,k).*sin(psi))/naz; 
else 

thetalc-2*sum(theta(:,k+1).*cos(psi))/naz; 
thetals-2*sum(theta(: ,k+l) .*sin(psi) ) /naz; 
end 

theta(:, k+1)-thetao+thetalc.*cos(psi)+thetals.*sin(psi); 

% *** calculation of new dthetadM *** 

theta-[theta theta(:,k+l)]; 

Mpsi = [Mpsi Mpsi (:, k+l)) ; 

theta(:. k+2)-theta(:,k)+0.25/57.3; 

Mpsi (: , k+2) -zeros (f^psi (: , k+1) ) ; 
k«k+2; 
tmcalc 
k-k-2; 

dthetadM-(theta(:, k+2)-theta(:,k*■./(Mpsi(:,k+2)-*$ssi(:,k)); 
k-k+1; 
end 


disp(' ') 

disp(' *** ADJUSTING COLLECTIVE ***') 

disp(' ') 

theta-theta(:,k); 
k—1; 

errorO-(TV 01)+1; 

while abs (errorO) > TV 01 
Tpsi-zeros(psi); 
thrcalc 

errorO«T-(mean(Tpsi)*b); 
if errorO < -TV 01, 

thetao-thetao-0.25*thetao*abs(l.25*errorO/T)*(l-mu); 
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elseif errorO > T*.01, 

thetao-thetao+O.25*thetao*abs(1.25*errorO/T)*(i-mu); 
•nd 

theta-thetao+thetalc.*cos(pai)♦chetala.*sin(pai); 
if k > 1, 

if aba(arrorO) > abs(errorl), 
clc 

diapt' ') 
diap(' ') 

disp ('This configuration will not trim*) 
dispt'Try a lower airspeed or a new design') 
disp(' ') 

disp('Program execution terminated') 
disp(* *) 

error('*** END OF PROGRAM ***') 
end 
end 

errorl-errorO; 
k-k+1; 
end 


% *** calculating drag moments **» 

DMpsi-zeros (psi) ; 
dmcalc 


% *** calculating rotor H force *** 

if Vinf < 16.9, 

Hrotor-0; 
dT- [dT ddT] ; 
dD«(dD ddD]; 
else 

dT* [dT ddT] ; 
dD- [dD ddD] ; 

for i«l:length(r)+l, 

Hlc(i)*2*sum(dT(:,i).*cos(psi))/naz; 

Hls(i)-2*sum(dD(:,i).*sin(psi))/naz; 
end 

Hrotor- (((b*cos (alphaT) /2) * (sum (His) - sin (betao) *sum(Hlc))) +Drotor) /2 ; 

end 


% *** calculating new rT *** 

rT2-(((mean(Mpsi(:,length(Mpsi(1,:))-1))/mean(Tpsi))/R)+rTl)/2; 


% *** check rotor drag and rT, retrim rotor if required *** 

while abs(Drotor-Hrotor) > 0.2*Hrotor | abs(rTl-rT2) > 0.015*rTl 

if abs(Drotor-Hrotor) > 0.2*Hrotor, 
disp(' ') 

disp(' *** ADJUSTING ROTOR DRAG ***') 

end 

Drotor-Hrotor; 
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if abs(rTl-rT2) > 0.015*zTl, 
disp(' ') 

dispC •** ADJUSTING MEAN THRUST LOCATION ***') 

•nd 


disp(' ') 

disp(' *** RETRIMMING ROTOR ***') 

dispC ') 
dTodT(:,1:nbe); 
dDmdD(:,1:nbe); 


% *** recalculating parameters *** 

alphaT*atan((Dftotal+Drotor)/(GW-Lftotal)); 
muoVinf*cos(alphaT)/Vtip; 

if Vinf >■ 16.9, 

T>(GW-Lftotal)/cos(alphaT); 

CT«T/(Adisk*rho*Vtip A 2); 

lamdaTo[i -2 w mu«sin (alphaT) 
mu A 2Msin (alphaT) A 2+l)-2*mu A 3*sin(alphaT) mu A 4*sin(alphaT) A 2-CT A 2/4]; 
lar -.T-max(real (roots (lamdaT))) ; 
vi= 'aT*Vtip-Vinf*sin(alphaT); 
viov*•ones(r); 
end 

B«l-(sqrt(2*CT)/b); 

Reff*B*R; 

RbaroReff-e; 

dro(Reff-grip)/nbe; 

rogrip:dr:Reff-dr;.ror+dr/2; 

RbarTorT2 *Rbar; 

betao«asin( (T/b*RbarT- (. c * (R-e) +e) *wblade) / ((.5* (R-e) -*-e) A 2*omega A 2*mblad 

e)) ; 


% *** trimming collective *** 

tOodock; 
k«l; 

errorOo(T*.02)+1; 

while abs(errorO) > T*.02 
Tpsiozeros(psi); 
thrcalc 

errorOoT-(mean(Tpsi)*b); 
if errorO < -T*.02, 

thetaoothetao-0.35*thetao*abs(1.5*errorO/T)*(1-mu); 
elseif errorO > T*.02, 

thetaoothetao+0.35*thetao*abs(1.5*errorO/T)*(1-mu); 
end 

thetaothetao+thetalc.*cos(psi)+thetals.*sin(psi); 
if k > 1, 

if abs(errorO) > abs(errorl), 
clc 

dispC ') 
dispC ') 
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(lisp('This configuration will not trim') 
disp('Try a lower airapaad or a now dosign') 
diap(' ') 

diapC Program oxocution terminated') 
diap(' ') 

error('*** END OF PROGRAM ***') 
and 
and 

errorl-error0; 
k-k+1; 
and 


% *** trimming cyclic *** 


k-1; 

errorO-(<(T/b)*rT2*(R-grip) )•.04)+1; 

while arrorO > ((T/b)*rT2*(R-grip))*.04 

time-etime(clock,tO); 
if time > 15, 

diep('still trimming ...') 
disp(' ') 
tO-clock; 
end 

Mpsi(:,k)-zeros(psi ); 
tmcalc 

theta-[theta theta(:,k)); 

Mpsi-[Mpsi Mpsi(:,k)j; 

% *** calculation of initial dthetadM *** 

if k < 2, 

theta(:,k+l)-theta(:,k)+0.25/57.3; 

Mpsi(:,k+1)-zeros(psi); 

k*k+l ; 

tmcalc 

k-k-l; 

dthetadM- (theta (: ,k+i) - theta (: ,k)) . / (Mpsi (: ,k*l) -I^psi (:, k)) ; 
end 

% *** calculation of M first harmonic parameters *** 

Mlc-2*sum(Mpsi(:,k).*cos(psi))/naz; 

Mls«2*sum(Mpsi(:,k).*sin(psi))/naz; 

% +** removal of first harmonic terms from J$psi *** 

Mpsi(:,k+l)-Mp8i(:,k)-Mlc.*cos(psi)-Mis.*ain(psi); 
delM-Mpsi(:,k+1)-Mpsi(:,k); 
errorO-max(delM)-min(delM); 
if k > 1, 

if errorO > errorl, 
clc 

disp(* ') 
disp(' ') 

disp('This configuration will not trim') 
disp('Try a lower airspeed or a new design') 
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dispC ') 

disp('Program execution terminated') 

disp(' ') 

error('*•* BUD OF PROGRAM **•') 
end 
end 

errorl-errorO; 

% *•* calculation of new theta **• 

delM-0.5*(1-mu)*delM; 

theta (:, k+l) -theta(: ,k) ♦ (dthetadM. *delM) ; 
if errorO <- ((T/b)*rT2*(R-grip))*.04, 
thetalc-2*sum(theta(:,k). *cos(psi))/naz; 
thetals-2*sum(theta(:,k).*sin(psi))/naz; 

else 

thetalc-2*sum(theta(:,k+l).*cos(psi))/nax; 
thetals-2+sum(theta(:,k+i).*sin(psi))/nax; 
end 

theta(:,k+l) -thetao+thetalc.*cos (psi) +thetals.»sin(psi) , 

% calculation of new dthetadM *** 

theta*[theta theta(:,k+1)]; 

Hpsi-[Mpsi Mpsi(:,k+1)]; 

theta(: , k+2)-theta(:,k)+0.25/57.3; 

Mpsi (: , k+2) -zeros (Mpsi (: , k+1) ) ; 
k«k+2; 
tmcalc 
k-k-2• 

dthetadM-(theta(:, k+2)-theta(:,k))./(Mpsi(:,k+2)-l$>si(:,k)); 
k-k+1; 
end 


% *** retrimming collective *** 

theta-theta(:,k); 
k—1; 

errorO-(T*.01)+l; 

while abs(errorO) > T*.01 
Tpsi-zeros(psi); 
thrcalc 

errorO-T-(mean(Tpsi)*b); 
if errorO < -T*.01, 

thetao-thetao-0.25*thetao*abs(1.25*errorO/T)*(1-mu); 
elseif errorO > T*.01, 

thetao-thetao+0.25*thetao*abs(1.25*errorO/T)*(1-mu); 
end 

theta-thetao+thetalc.*cos(psi)+thetals.*sin(psi); 
if k > 1, 

if abs(errorO) > abs(errorl), 
clc 

disp(' ') 
disp(' ') 

disp('This configuration will not trim') 
disp('Try a lower airspeed or a new design') 
disp(* ') 

disp('Program execution terminated') 
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diap (' ') 

error('*** END OF PROGRAM **•*) 
•nd 
•nd 

•rrorl-errorO; 
k-k+1; 

•nd 


% *** recalculating rotor H force *** 

if Vinf < 16.9, 

Hrotor-0; 
dT- IdT ddT]; 
dO- IdD ddD] ; 

•lee 

dT- [<JT ddT] ; 
dO-[dO ddD]; 

for i-l:length(r)+1, 

Hlc(i)-2*sum(dT(:,i). *cos(psi))/naz; 

Hls(i)*2*sum(dD(:,i),*sin(psi))/naz; 
end 

Hrotor- (((b*cos (alphaT) /2) * (sum (Hla) - sin (betao) *Bum(Hlc))) -*-Drotor) /2 ; 
end 


% *** recalculating rT *** 

rTl*rT2; 

rT2* ( ( (mean (l^pBi {:, length U^>si (1, :)) -1)) /mean (Tpsi)) /R) +rTl) /2; 
end 


% *** recalculating drag moments *** 

dT*dT(:,1:nbe); 
do-dO (:, l ,-nbe) ; 

Dl^psi-zeros (psi) ; 

dmcalc 

dT- [dT ddT] ; 

dD- [dO ddD]; 

clc 

disp(' ') 
diap(' ') 

dispC *** ROTOR TRIMMED ***') 

pause(3) 
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c 


CRUISE.* 


» CRUISE.M 
% CALLED BY STAB.M 

% Computes the stability derivatives in cruise flight. 

% calls the following subroutines 

» 

% CMRGRP 
% CTRGRP 
% CFUSEGRP 
% TRIM 

% DCTPLOTS or DCTMATS 

% 

% evaluate basic stability derivatives 

% 

% MAIN ROTOR 

% dctsigdmu dctsigdtheto dctsigdlamp 

% dchsigdmu dchsigdtheto dchsigdlamp 

% dcqsigdmu dcqsigdtheto dcqsigdlamp 

% daldmu daldtheto daldlamp 

% dbldmu dbldtheto dbldlamp 

% 

% TAIL ROTOR 

% dctsigdlampt dctsigdmut dctsigdthetot 

% 

vctsig-ones(1,3); 
vcqsigsones(1,3); 
vchsig-ones(1,3); 
valtpp-ones (1,3) ; 
vAl«ones(1,3) ; 
vBl»ones(l,3); 
vals*ones(1,3); 
vbls«ones(l,3) ; 
vmusones(1,3) ; 
vtheta7-ones(1,3); 
vao»ones(1,3); 
wl*ones (1,3) ; 
vlamp=ones(1,3); 
vthetao*ones(1,3); 
vctsigt*ones(1,3); 
vmutsones(1,3); 
vlampt=ones(1,3); 
vthetaotsones(1,3); 

% 

thetal*-twist; 
eval(['load ',filenamel)) 

Vinf*Vinf*1.0005; % .05% higher 

disp(' ') 

disp(' *** FIRST OF THREE TRIM ITERATIONS ***') 

disp(' ') 
pause(3) 
trim 

Qrotor-mean(DMpsi)*b; 
solidity«b*cblade/(pi*R); 

CQ-Qrotor/(Adisk*rho*Vtip“2*R); 

CH*Hrotor/(Adisk*rho*Vtip A 2); 
vctsig(1,3)■CT/solidity; 
vcqsig(1,3)»CQ/solidity; 
vchsig(l,3)«CH/solidity; 
valtpp(1,3)--alphaT; 
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vAl(1,3)--thetalc; 

vBl(1,3)--thetal8; 

vmu(l,3) -mu; 

thetaiave-thetao; 

vthetao(1,3)-thetaaave-.7 *thecal; 

vthata7(1,3)-vthetao(1,3)♦.7*thetal; 

vao(l,3)-betao; 

wl (1,3) -inaan(vi) ; 

vlamp (1,3) -vmu(l,3) *valtpp(l,3) - wl (1,3) /omega/R; 

lockno-rho*a*cblade*R A 4/Ib; 

numl-vthetao(l,3)*(8/3+32/45*vmu(1,3 )“2 /pi) ; 

num2-thetal*(2+vmu(l,3) A 4/12); 

num3-vlampd, 3)*(2-vmu(1,3)*vmu(1,3)/2); 

denl-l+3/2*vmu(l,3) *vmu(l,3) -5*vmu(l,3) *4/24; 

vala(l,3)-vmu(l,3)*(numl+num2+num3) /denl-vBl (1,3) ; 

vbls(1,3)-vals(1,3)*l2/lockno*e/R/(l+e/R/3); 

eqsig-CQ/solidity; 

ohm-omega; 

eigma-solidity; 

lv*lvd-xcg; 

g«32.2; 

A-pi*R*R; 

» 

% tail rotor 

% 

if notar--0 
Abt*bt*Rt*cot; 

At*pi*Rt*Rt; 

mut-Vinf/ohmt/Rt; 

sigmat-Abt /At; 

lt-ltd-xcg; 

ht-htd-rcg; 

yt*ytd-ycg; 

locknot«rho*at*cot*Rt*4/Ibt; 

Tt=<cqsig*sigma*rho*A*(ohm*R) A 2*R-Lvert*lv)/It; 
ctsigt-Tt/(rho*At*(ohmt*Rt)*2); 
lampt--ctsigt*sigma/2/mu; 

aot*2/3*locknot*ctsigt/at-3/2*g*Rt A 2/(ohmt*Rt)*2; 

alst-(-2*((4/3*mut*aot+ctBigt*sigmat/2/mut)/(2+mut*mut))*(l-(4*mut/(2+3* 
mut*mut)) A 2)*... 


tan(delta3))/(((2-mut*mut)/(2+3*mut*mut))+(1-(4*mut/(2+3*mut*mut)) A 2)*ta 
n(delta3) A 2); 

bl8t-2*((4/3*mu*aot*ctsigt*sigma/2/mu)/(2+mu*mu))+alst*tan(delta3); 

thetaot-((4*ctsigt/at)-(.5+mut*mut/2)*thetalt+mut*blst*tan<delta3)-lampt 

)/... 

(2/3+mut*mut)-aot*tan(delta3); 
theta75t-thetaot+.75*thetalt; 
vlt-sqrt(Tt/2/rho/At); 
vctsigt(1,3)-ct8igt; 
vmut (1,3)-mut; 
vthetaot(1,3)-thetaot; 
vlampt (1,3) -lampt; 
end 
% 

aave stabtemp vctsig vcqsig vchsig valtpp vAl vBl vmu vthetao thetal ... 
vtheta7 vao wl vlamp lockno vals vbls rho filenamel vctaigt vmut ... 
vthetaot vlampt g 
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clear 

load atabtanp 

aval(['load ',filanamel]) 

load stabtemp 

Vinf*Vinf*.9995; % .05% lower 

V.Vinf; 
disp(' ') 

disp(' *** SECOND OF THREE TRIM ITERATIONS ***') 

di»p(' ') 
pause(3) 
trim 

Qrotor*mean(DMpsi)*b; 
solidity*b*cblade/(pi*R); 

CQ*Qrotor/(Adisk*rho*Vtip A 2*R); 

CH*Hrotor/(Adisk*rho*Vtip A 2); 

vetsig(1,1)*CT/solidity; 

vcqsigd, 1) *CQ/solidity; 

vchsig(1,1)*CH/solidity; 

valtpp(1,1)*-alphaT; 

vAl (1,1)*- thetalc,- 

vBl(1,1)*-thetals; 

vmu(1,1)*mu; 

thetasave*thetao; 

vthetao(1,1)*thetasave-.7*thetal; 

vtheta7 (i,i) *vthetao(l, 1) +.7*thetal ; 

vao(1,1)*betao; 

wl (1,1) *mean(vi) ; 

vlamp(1,1)*vmu(1,1)*valtpp(1,1)-wl(1,1)/omega/R; 
lockno*rho*a*cblade*R A 4/Ib; 

numl*vthetao(1,1)*(8/3+32/45*vmu(1,1) A 3/pi); 

num2*thetal*(2+vmu(l,1) A 4/12); 

num3*vlamp(1,1)*(2-vmu(l,1)*vmu(l,1)/2); 

denl*l+3/2*vmu(l, 1) *vmu(1,1) -5*vmu(1,1) *4/24 ,- 

vale (1,1) »vmu (1,1) * (numl+num2+num3) /denl-vBl (1.1) ; 

vbls(1,1)*vals(1,1)*12/lockno*e/R/(l+e/R/3); 

cqsig*CQ/solidity; 

ohm*omega; 

sigma*solidity; 

lv*lvd-xcg; 

A*pi*R*R; 

% 

%tail rotor 

% 

if notar**0 
Abt=bt*Rt*cot ; 

At*pi*Rt*Rt; 
mut*Vinf/ohmt/Rt; 
sigmat*Abt/At; 
lt*ltd-xcg; 
ht-htd-zcg; 
yt*ytd-ycg; 

locknot*rho*at*cot*Rt A 4/Ibt; 

Tt*(cqsig*sigma*rho*A*(ohm*R) A 2*R-Lvert*lv)/It; 
ctsigt*Tt/(rho*At*(ohmt*Rt) A 2); 
lampt*-ctsigt*sigxna/2/mu; 

aot*2/3*locknot*ctsigt/at-3/2*g*Rt A 2/(ohmt*Rt) A 2; 

alst*(-2* ((4/3*mut*aot+ctsigt*sigmat/2/mut) / (2+mut*mut))*(1- (4*mut/ (2+3* 
mut*mut)) A 2)*... 
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tan(delta3)) / (((2-mut*mut) / (2+3*mut*mut) ) ♦ (1- (4*mut/ (2+3*mut*mut)) A 2)*ta 
n(delta3)*2); 

blst«2*((4/3*mu*aot*ctsigt*sigma/2/mu)/(2+mu*mu))+alst*tan(d«lta3); 

thetaot*((4*ctsigt/at)-(.5+mut*mut/2)*thetalt-Hnut*blst*tan(delta3)-lampt 
)/••- 

(2/3+mut*mut)-aot*tan(delca3); 
theta75tathetaot+.75*thetalt; 
vltasqrt(Tt/2/rho/At); 
vctaigt(1,1)actsigt; 
vmut(1,1)smut; 
vthetaot(1,1)«thetaot; 
vlampt(1,1)alampt; 

•nd 

% 

save stabtemp vctsig vcqsig vchsig valtpp vAl vBl vmu vthetao thetal ... 
vtheta7 vao wl vlamp lockno vals vbls rho filenamel vctsigt vmut ... 
vthetaot vlampt g 
clear 

load stabtemp 
eval(['load filenamel]) 
load stabtemp 
disp(' ') 

dispC *** THIRD OF THREE TRIM ITERATIONS ***') 

disp(* ') 
pause(3) 
trim 

Qrotoramean (DMpsi) *b; 
solidity=b*cblade/(pi*R); 

CQ=Qrotor/ (Adisk*rho*Vtip A 2*R) ; 

CHaHrotor/(Adisk*rho*Vtip A 2); 
disp(' ') 

disp(' *** PERTURBATIONS COMPLETE ***') 

disp(' ') 

disp(' *** EVALUATING STABILITY DERIVATIVES ***') 

disp(' ') 

vctsig(l,2)aCT/solidity; 

vcqsig(1,2)aCQ/solidity; 

vchsig(1,2)aCH/solidity; 

valtpp(1,2)=-alphaT,- 

vAl(1,2)a-thetalc; 

vBl(l,2)a-thetals; 

vmu(1,2)*mu; 

the tasave a the tao ; 

vthetao(1,2)athetasave-.7*thetal; 

vtheta7(1,2)avthetao(l,2)♦,7*thetal; 

vao(1,2)abetao; 

wl (l,2)amean(vi) ; 

vlamp (l, 2) avmu (1,2) *valtpp (1,2) - wl (1,2) /omega/R; 

V-Vinf; 

locknoarho*a*cblade*R A 4 /lb ,- 

numlavthetao(1,2)*(8/3+32/45*vmu(1,2)*3/pi); 
num2athetal*<2+vmu(1,2) A 4/12); 
num3avlamp(l,2)*(2-vmu(1,2)*vmu(1,2)/2); 
denlai+3/2*vmu(l,2)*vmu(1,2)-5*vmu(1,2) A 4/24; 
vals(1,2)avmu(l,2)*(numl+num2+num3)/denl-vBl(1,2); 
vbls(1,2)avals(1,2)*12/lockno*e/R/(l+e/R/3); 
cqsig*CQ/solidity; 
ohmaomega; 
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sigma-solidity; 

lvalvd-xcg; 

A»pi*R*R; 

% 

% tail rotor 

% 

if notarssO 
Abt-bt*Rt*cot; 

At»pi*Rt*Rt; 

mut «Vinf /ohmt /Rt ; 

sigmataAbt/At; 

lt-ltd-xcg; 

htahtd-zcg; 

yt«ytd-ycg; 

1ocknot«rho*at*cot*Rt A 4/Ibt; 

Tt« (cqsig*sigma*rho*A* (ohm*R) A 2*R-Lvert*lv) /It; 
ctsigt«Tt/(rho*At*(ohmt*Rt)“2); 

1 anqpt a-ctsigt* sigma/2/mu; 

aota2/3*locknot*ctsigt/at-3/2*g*Rt A 2/(ohmt*Rt) A 2; 

alst* (-2* ((4/3*mut*aot+ct6igt*sigmat/2/mut) / (2+mut*mut)) * (1- (4*mut/ (2+3* 
mut*mut))*2)*... 

t«m (delta3)) / ( ((2-mut*mut)/(2+3*mut*mut)) + (1-(4*mut/(2+3*mut*mut)) A 2)*ta 
n(delta3 ) a 2); 

blst*2*((4/3*mu*aot*ctsigt*sigma/2/mu)/(2+mu*mu))+alst*tan(delta3); 

thetaot*((4*ctsigt/at)-(.5+mut*mut/2)*thetalt+mut*blst*tan(delta3)-lampt 

)/... 

(2/3+mut*mut)-aot*tan(delta3); 
theta75t»thetaot+.75*thetalt; 
vlt=sqrt(Tt/2/rho/At); 
vetsigt(1,2)actsigt; 
vmut(1,2)=mut; 
vthetaot(1,2)athetaot; 
vlampt(1,2)*lampt; 
end 
% 

save stabtemp vetsig vcqsig vchsig valtpp vAl vBl vmu vthetao thetal ... 
vtheta7 vao wl vlantp lockno vals vbls T rho filenamel 

% 

% Solving for unknowns 

% 

DCTPLOTS 

% 

dcysigdbxinv(T*sin(.1)); 

% 

if notaraaO 

dctplott 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%»%%*%%%%%%%% 

% CONFIGURATION CALCULATION 

% 

% converting data from crank 
theta7avthetao(l,2) 
thetaoatheta7-.7*thetal; 
ctsigavctsig(1,2); 
chsig«vchsig(1,2); 
cqsigavcqsig(1,2); 
altppavaltpp(1,2); 
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Al.vAl(1,2); 

Bl-vBl(l,2) ; 
aOxvaoU, 2) ; 
vlxwl (1,2) ; 

VxVinf ; 

t Main rotor 

hmxhmd-scg; 

ym*ymd-ycg; 

lntxlmd-xcg; 

c-cblade; 

Ab»c*R*4; 
lampxvlamp(l,2) ; 
mu-vmu(1,2); 
ctxCT; 

fxAfV; 

thetatx4/a*ctsig+sqrt(sigma*ctsig/2); 

lockno«rho*a*c*R A 4/Ib; 

theta75«thetao+.75*thetal; 

alsovals(l,2); 

blsxvbls(1,2); 

m»GW/g; 

Zc«Xxx*Xzz; 
lhslhd-xcg; 
hh»hhd-zcg; 
yh«yhd-ycg; 
hvxhvd-zcg; 
yvxyvd-ycg; 

delvmaxxpi/4; test max verticle fin defl before stall 
if maxr>>0 
ddelvddelp*0; 
else 

ddelvddelpxdelvmax/maxr; % verticle fin deflection/maxr 
end 

if notarxxi; 

% NOTAR 

vlmaxxi,2*vi; % downwash at the boom slots (NOTAR) 

ltnxlttnd-xcg; 

lnxltnd-xcg; 

ynxytnd-xcg; 

hnxhtnd-zcg; 

htnxhn; 

Tt«(cqsig*sigma*rho*A*(ohm*R) A 2*R-Lvert*lv)/ltn; 

cmunx.55; % a typical value 

bn«.5*R; 

hslotx.028*dian/2; 
arnxbn/dian; 

At«pi*dian*dian/4 ; 

Rtxdian/2; 
end 

if Swing<.1 
wingxO ; 
else 

AwxSwing; 

lwxlwd-xcg; 

hwxhwd-zcg,- 

ywxywd-ycg; 

alpwxCLwing/aw; 

deliwxewing; 

ZwxLwing; 
bwxbwing; 
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cdowaCDowing; 
wingsl; 

•nd 

Av.Svert; 
bvabvert; 

Ah-Shoriz; 
bhabhoriz; 
cdov*CDovert; 
cdohaCDohoriz; 
delivaewing; 
delih*ewing; 
alphaCLhoriz/ah; 

ZhaLhoris; 

YvaLvert; 

if Yv<.1;Yva. 1; e id; 

deldval.72*(1-b ' (20+Rt))*(2*Yv/(2*Rt+bv)); 

% trim conditions 
thoa-T*(im+als)/GW; 
phoa-(T*bls+Tt)/GW; 
wo*0; 

VOaO ; 
uoaVinf; 
games0; 

%%%%%%%%%%%%%%%%%«%%«%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% STABILITY CALCULATION 

» 

cmrgrp 

ctrgrp 

cbodygrp 

% 

% computation of A,B,C,D matrices 

% 

% CHECK THE WING CALCULATIONS 

% 

All* (dxdxdotm+dxdxdoth+dxdxdotv+dxdxdotf+dxdxdotw)/m; 

A12* (dxdzdotm+dxdzdoth+dxdzdotf+dxdzdotw)/m; 

A13= (dxdqm)/m-wo; 

A14 = -g*cos(tho); 

A15* (dxdydotm+dxdydotv)/m; 

A16* dxdpm/m; 

A17* 0; 

A18 = 0;%(dxdrw)/m+vo; recheck this one 

% 

A21* (dzdxdotm+dzdxdoth+dzdxdotf+dzdxdotw)/m; 

A22* (dzdzdotm+dzdzdoth+dzdzdotf+dzdzdotw)/m; 

A23* (dzdqh+dzdqw)/m+uo; 

A24* -g*cos(pho)*sin(tho); 

A25* 0;% (dzdydotw)/m; 

A26* 0;% (dzdpw)/m-vo; 

A27* -g*sin(pho)*cos(tho); 

A28a dzdrm/m;% (dzdrm+dzdrw)/m; % recheck this one 

% 

A31* (dmdxdotm+dmdxdoth+dmdxdotf+dmdxdotw)/Iyy; 

A32a (dmdzdotm+dmdzdoth+dmdzdotf+dmdzdotw)/Iyy; 

A3 3■ (dmdqm+dmdqh+dmdqw)/Iyy; 

A34- 0; 

A35a dmdydotm/Iyy; 

A36- (dmdpm)/Iyy; 

A37- 0; 
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A38« 0;% dmdrf/Iyy; 

t 

A41- 0; 

A42- 0; 

A43« cos(pho); 

A44- 0; 

A45- 0; 

A46- 0; 

A47« 0; 

A4B« -sin(pho); 

% 

A51<> (dydxdotm+dydxdott+dydxdotn+dydxdotv) /n; 

A52- (dydzdotm+dydzdotn)/m; 

A53- dydqm/m; 

A54b -g*sin(pho)*sin(tho); 

ASSb (dydydotm+dydydott+dydydotn+dydydotv+dydydotf) /■; 

A56« (dydpm+dydpt+dydpn+dydpv)/m+wo; 

A57* g*cos(pho)*COS(tho); 

A58« (dydrt+dydrn+dydrv)/m-uo; 

% 

A61» (Izz*(drdxdotm+drdxdott+drdxdotn+drdxdotv)... 

♦ Ixz*(dndxdotm+dndxdott+dndxdotn+dndxdotv))/Zc; 

A62» (Izz*(drdzdotnw-drdzdotn)+Ixz*(dndzdotm+dndzdotn))/Ic; 

A63« (Izz* (drdcjm) +Ixz* (0)) /Ic; % recheck this one 
A64« 0; 

A65» (Izz*(drdydotm+drdydott+drdydotn+drdydotv+drdydotf)... 

♦Ixz*(dndydott+dndydotn+dndydotv+dndydotf))/Ic; 

A66* (Izz*( drdpm+drdpt+drdpn+drdpv+drdpw)tlxi*(dndpt+dndpn+dndpv))/Ic; 
A67* 0; 

A68* (Izz*(drdrt+drdrn+drdrv+drdrw)+lxz*(dndrm+dndrt+dndm+dndrv))/Ic; 

% 

A71- 0; 

A72« 0; 

A73« sin(pho)*tan(tho); 

A74* 0; 

A75* 0; 

A76» 1; 

A77- 0; 

A78* cos(pho)*tan(tho); 

% 

A81* (Ixz*(drdxdotm+drdxdott+drdxdotn+drdxdotv)... 

+Ixx*(dndxdotm+dndxdott+dndxdotn+dndxdotv))/Ic; 

A82* (Ixz* (drdzdotm+drdzdotn) +Ixx* (dndzdotm+dndzdotn) ) /Ic; 

A83« (Ixz*(drdqm)+Ixx*(0))/Ic;% recheck 
A84* 0; 

A85» (Ixz*(drdydotm+drdydott+drdydotn+drdydotv+drdydotf)... 

♦Ixx*(dndydott+dndydotn+dndydotv+dndydotf))/Ic; 

A86» (Ixz*(drdpm+drdpt+drdpn+drdpv+drdpw)+Ixx*(dndpt+dndpn+dndpv))/ic; 
A87* 0; 

A88*(Ixz*(drdrt+drdrv+drdrw)+Ixx*(dndrt+dndrv))/Ic; 

» 

% longitudinal plant augmented is X«[u w q theta]' 

Flonaug* [All A12 A13 A14 ;A21 A22 A23 A24;A31 A32 A33 A34 ; 

A41 A42 A43 A44] ; 

Plonaug-poly(Flonaug); 

Rlonaug-roots(Plonaug); 

% lateral plant augmented is X-[v p phi r psi]' 

Flataug* (A55 A56 A57 A58 0;A6S A66 A67 A68 0; 

A7S A78 A77 A78 0;A85 A86 A87 A88 0;0 0 0 1 0] ; 

Plataug-poly(Flataug); 





Rlataug«roots(Plataug); 

% coupled plant 

Amat» [All A12 A13 A14 A15 A16 A17 A18 0; 

A21 A22 A23 A24 A25 A26 A27 A28 0; 

A31 A32 A33 A34 A35 A36 A37 A38 0; 

A41 A42 A43 A44 A45 A46 A47 A48 0; 

A51 A52 A53 A54 A55 A56 AS7 A58 0; 

A61 A62 A63 A64 A65 A66 A67 A68 0; 

A71 A72 A73 A74 A75 A76 A77 A78 0; 

A81 A82 A83 A84 ASS A86 A87 A88 0; 

00000001 0 ); 

Pcoup>poly(Amat); 

Rcoup-roots(Pcoup); 

» 

Bllo dxdbim* dblmdde1e/m; 

B12- dxdthetotn*dthetomddelc/m; 

B13« dxdalm*dalmddela/m; 

B14* 0; 

% 

B21* dzdblm*dblmddele/m; 

B22« dzdthetom*dthetomddelc/m; 

B23- 0; 

B24« 0; 

% 

B31» dmdblm*dblmddele/Iyy; 

B32« dmdthetom*dthetomddelc/Iyy; 

B33* dmdalm*dalmddela/Iyy; 

B34- 0; 

% 

B41« 0; 

B42> 0; 

B43* 0; 

B44- 0; 

% 

B51« dydblm*dblmddele/m; 

B52» dydthetom*dthetoir.ddelc/m; 

B53= dydalm*dalmddela/m; 

B54= (dydphin*dphinddelp+dyddelv*ddelvddelp+dydthetot*dthetotddelp)/m; 

% 

B61* Izz*drdblm*dblmddele/Ic; 

B62« (Izz*drdthetom*dthetomddelc+Ixz*dndthetom*dthetomddelc)/Ic; 

B63= Izz*drdalm*dalmddela/Ic; 

B 6 4 

(Izz*(drdphin*dphinddelp+drddelv*ddelvddelp+drdthetot*dthetotddelp)+... 

Ixz*(dndphin*dphinddelp+dnddelv*ddelvddelp+dndthetot*dthetotddelp))/Ic; 

% 

B71« 0; 

B72- 0; 

B73- 0; 

B74« 0; 

% 

B81« Ixz*drdblm*dblmddele/Ic; 

B82« (Ixz*drdthetom*dthetomddelc+Ixx*dndthetccn*dthetofnddelc) /Ic; 

B83o Ixz*drdalm*dalmddela/Ic; 

B 8 4 

(Ixz*(drdphin*dphinddelp+drddelv*ddelvddelp+drdthetot*dthetotddelp)+... 

Ixx*(dndphin*dphinddelp+dnddelv*ddelvddelp+dndthetot*dthetotddelp))/Ic; 

% 








Glonaug*[Bll B12 B13 B14; 

B21 B22 B23 B24 ; 

B31 B32 B33 B34; 

B41 B42 B43 B44]; 

Glataug-[BS1 B52 B53 B54; 

B61 B62 B63 B64; 

B71 B72 B73 B74; 

B81 B82 B83 B84; 

0 0 0 0 ]; 

% coupled input matrix 

Bmat«[Bll B12 B13 B14.B21 B22 B23 B24; 

B31 B32 B33 B34;B41 B42 B43 B44; 

B51 B52 B53 B54.B61 B62 B63 B64; 

B71 B72 B73 B74.B81 B82 B83 B84.0 000]; 

% 

xcouple«12/lockno*e/R/(l*e/3/R); 

% designed damping 
desdmdq«dmd<pn+dmdqh+dmdqv ; 
desdrdp*drdpm+drdpt+drdpv+ drdpv+drdpn ; 
deBdndr-dndrm+dndrt+dndrv+dndrn; 

% cooper harper pilot rating 
prpitch-desdmdq/Iyy; 
prrol1*de sdrdp/Ixx; 
pryaw«desdndr/Izz; 

% control power 
cppitch*B3l*lyy; 
cproll«B63*Ixx; 
cpyaw«B84 *Izz; 
cpipitch*B31; 
cpiroll»B63; 
cpiyaw.B84; 

% 

theta0-theta7; 
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D 


DCTPLOTS.M 


% DCTPLOTS.M 

% 

% plots ct/sigma, ch/sigma, cq/sigma, als and bis 
% VERSES 

% mu, thstao and lamp 

% fits a polynomial to the curve, takes a derivative of the curves 
t then evaluates the appropriate parameters for CRUISE.M 

I 

t d(ct/sigma), d(ch/sigma), d(cq/sigma), dais and dbls 
» VERSES 

% dmu, dthetao and dlamp 

% uses the pertubation points from CRANK.M about nominal 

% 

% polyfitting data to obtain plots for n-8 
n«2; 

ectsigmu-polyfit(vmu,vctsig,n); 
ectsigtheto-polyfit(vthetao,vctsig,n),- 
ectsiglamp.polyfit(vlamp,vctsig,n); 

% 

echsigmu-polyfit(vmu,vchsig,n); 
echsigtheto*polyf it (vthetao, vchsig,n) ,- 
echsiglamp*polyfit(vlamp,vchsig,n); 

% 

ecqsigmu*polyfit(vmu,vcqsig,n); 
ecqsigtheto*polyfit(vthetao,vcqsig.n),- 
ecqsiglamp*polyfit(vlamp,vcqsig.n); 

% 

ealsmuspolyfit(vmu,vals.n); 
ealstheto=polyfit(vthetao,vals.n); 
ealslamp*polyfit(vlamp,vals.n); 

% 

eblsmu^polyfit(vmu,vbls.n); 
eblstheto=polyfit(vthetao,vbls,n); 
eblslamp-polyfit(vlamp,vbls.n); 

% 

%taking derivatives 

% 

edctsigdmu=deriv2(ectsigmu); 
edctsigdtheto*deriv2(ectsigtheto); 
edctsigdlamp*deriv2 (ectsiglamp) ; 

% 

edchsigdmu*deriv2(echsigmu); 
edchsigdtheto=deriv2(echsigtheto); 
edchsigdlamp=deriv2(echsiglamp); 

% 

edcqsigdmu*deriv2(ecqsigmu); 
edcqsigdtheto*deriv2(ecqsigtheto); 
edcqsigdlamp*deriv2(ecqsiglamp); 

% 

edalsdir.^-deriv2 (ealsmu) ; 
edalsdtheto*deriv2(ealstheto); 
edalsdlamp=deriv2 (ealslamp) ; 

% 

edblsdmu«deriv2(eblsmu); 
edblsdtheto«deriv2(eblstheto); 
edblsdlamp«deriv2 (eblslamp) ; 

% 

% evaluating derivatives 


137 



% 

dct*igdmu»polyv«l(adctaigdmu,vmu(l,2)); 
dccaigdcheto»polyval(•dct»igdth«to,vth«taotl / 2)); 
dctaigdlampapolyval (adctsigdlanp, vlainp (l, 2)) ; 

» 

dchaigdnu-polyval(edchaigdmu,vmu(1,2)); 
dchaigdthato-polyval(adchaigdthaco,vthatao(1,2)); 
dchsigdlanpapolyval(adchaigdlamp,vlamp(1,2)); 

» 

dcqaigdmu»polyval(adcqaigdmu,vmu(1.2)); 
dcqaigdthato-polyval(adcqaigdthato,vthatao(1,2)); 
dcqsigdlanp«polyval (adcqnigdlamp, vlainp (1,2)) ; 

% 

daldmu«polyval(adalsdmu, vmu(1,2)); 
daldthato-polyval(adaladthato,vthatao(1,2)); 
daldlanp-polyval (adaiadlamp,vlamp(l,2)) ; 

% 

dbldmu«polyval(adbladmu,vmu(1,2)); 
dbldthetoapolyval(edblsdtheto,vthatao(1,2)); 
dbldl&mp-polyval(adbladlamp,vlamp(1,2)); 

% 
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S. DCTH»OTT.It 


» DCTPLOTT.M 

% 

% plot* ct/sigmat 
« VERSES 

% mut, thetaot and lampt 

% fits a polynomial to the curve, takes a derivative of the curves 
t then evaluates the appropriate parameters for CRUISE.M 

% d(ct /sigma) 

% VERSES 

% dmut, dthetaot and dlampt 

% uses the pertubation points from JANRAD about nominal 

% 

% polyfitting data to obtain plots for n>2 
n>2; 

ectsigmutspolyfit(vmut,vctsigt.n); 
ectsigthetotapolyfit(vthetaot,vctsigt,n); 
ectsiglampt«polyfit(vlampt,vctsigt.n); 

% 

ttaking derivatives 

% 

edctsigdmut»deriv2(ectsigmut); 
edctsigdthetot»deriv2(ectsigthetot); 
edctsigdlampt*deriv2(ectsiglampt); 

% 

% evaluating derivatives 

» 

dctsigdmut-polyval(edctsigdmu,vmu(1,2)); 
dctsigdthetot-polyval(edctsigdtheto,vthetao(1,2)); 
dctsigdlampt=polyval(edctsigdlamp,vlamp(1,2)); 
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r. ama&r.M 

t cmrqrp.M 
% CALLED BY CRUISE.M 

% Computes the basic MR derivatives in cruise flight 

% Computes the stability derivatives of the main rotor in cruise flight. 
% Uses data loaded in the workspace by JANRAD.M and STAB.M 

% 

% Compute the basic mainrotor derivatives. 

I 

dmudxdotai/ohm/R; 

dlampdxdotai/ohm/R*(altpp-sigma/2/mu*(dctsigdmu-ctsig/mu)); 

dlampdzdotal/(ohm*R*(l+dctsigdlamp*sigma/2/mu)); 

dbetdydotai/V; 

dchsigdaactsig+a/8*lamp; 

dchsigdbadchsigda; 

daldqa-16/<lockno*ohm*(1-e/R) A 2* (l-mu*mu/2))-... 

12*e/R/(lockno* ohm*(1-e/R) A 3*(l-mu A 4/4)); 

daldp-1/(ohm*(l-mu A 2/2))-192*e/R/(lockno A 2*ohm*(1-e/R) A 5*(1-mu A 4/4)); 
daldaal2*e/R*(l+mu A 2/2)/(lockno*(1-e/R) A 3*(l-mu A 4/4)); 
daldba-1*(l+3/2*mu A 2)/(l-mu A 2/2); 

dbldqa-inv(ohm*(l+mu A 2/2))+192*e/R/(lockno A 2*ohm*(1-e/R) A 5*(l-mu A 4/4)); 
dbldpa-16/(lockno*ohm*(1-e/R) A 2*(l+mu*mu/2))-... 

12*e/R/(lockno*ohm*(1-e/R) A 3*(1-mu A 4/4)); 
dbldaal; 

dbldbal2*e/R*(l+3/2*mu A 2)/(lockno*ll-e/R) A 3*(l-mu A 4/4)); 
dmdalsma3/4*e/R*Ab*rho*R*(ohm*R) A 2*a/lockno; 
drdblsmadmdalsm; 

% 

% Compute the mainrotor stability derivatives. 

% 

dxdxdotma-rho*Ab*(ohm*R) A 2*((dchsigdmu+dchsigda*daldmu+(als+im)*dctsigdm 

u) . . . 

*dmudxdot+((dchsigdlamp+dchsigda*daldlanp+(als+im)*dctsigdlanp)*dl«uipdxd 
ot)) ; 

dxdydotmarho*Ab*(ohm*R) A 2*ctsig*(Al-bls)*dbetdydot; 

dxdzdotma-rho*Ab*(ohm*R) A 2*(dchsigdlamp+ctsig*daldlamp+(als+im)*dctsigdl 
amp)... 

*dlampdzdot; 

dxdqma-rho*Ab*(ohm*R) A 2*dchsigda*daldq-dxdxdotm*hm; 
dxdpma-rho*Ab*(ohm*R) A 2*dchsigda*daldp-dxdydotm*hm; 
dxdthetoma-rho*Ab*(ohm*R) A 2*(dchsigdtheto+dchsigda*daldtheto+... 

(als+im)*dctsigdtheto); 

dxdalma-rho*Ab*(ohm*R) A 2*dchsigda*dalda; 
dxdblma-rho*Ab*(ohm*R) A 2*dchsigda*daldb; 

dydxdotmarho*Ab*(ohm*R) A 2*(dcysigdb*dbldmu+bls*dctsigdmu)*dmudxdot; 
dydydotma-rho*Ab*(ohm*R) A 2*(chsig+ctsig*(Bl+als))*dbetdydot; 
dydzdotmarho*Ab*(ohm*R) A 2*dbldlanp*dcysigdb*dlanpdzdot; 
dydqma rho*Ab*(ohm*R) A 2*dcysigdb*dbldq+dydxdotm*hm; 
dydpmarho*Ab*(ohm*R) A 2*dcysigdb*dbldp+dydydotm*hm; 
dydthetcmarho*Ab*(ohm*R) A 2*(dcysigdb*dbldtheto+bls*dctsigdtheto); 
dydalm«rho*Ab*(ohm*R) A 2*dcysigdb*dblda; 
dydblmarho*Ab*(ohm*R) A 2*dcysigdb*dbldb; 

dzdxdotma-rho*Ab* (ohm*R) A 2* (dctsigdmu*dmudxdot+dctsigdlazip*dlaiipdxdot) ; 

dzdzdotma-rho*Ab*(ohm*R) A 2*dctsigdlamp*dlaflpdzdot; 

dzdrma2/ohm*rho*Ab*(ohm*R) A 2*ctsig; 

dzdthetoma-rho*Ab*(ohm*R) A 2*dctsigdtheto; 

dzdblma-rho*Ab*(ohm*R) A 2*dctsigdlamp*inv(daldlanp)*daldb; 

drdxdotmadrdblsm*dbldmu*dmudxdot+dydxdotm*hm+dzdxdotm*ym; 

drdydotma-drdblsm*(Bl+als)*dbetdydot+dydydotm*hm; 
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drdsdotn>drdblsa*dbldlainp*dlainpdzdot+dydzdotin*hm+dzdzdotm*ym; 

drdqm»drdbl8m*dbldq+dydqm*hxn; 

drdpni»drdbl8m*dbldp+dydp!n*hm ; 

drdthe tom-dr dblBin*dbldtheto+dydthetom«hm+dzdthetoni*yin; 
drdalm-drdbl am* dblda+dydalm*hm; 
drdblm-drdblsm*dbldb+dydblm*hm; 

dmdxdotm-dmdalsm*(daldmu*dmudxdot+daldlamp*dlainpdxdot)+dxdxdotm*hm+dzdxd 

dmdydotm»-dxdydotm*hm+dn»dalBm* (Al-ble) *dbetdydot; 
dmdzdotm-dmdal8m*daldlainp*dlajnpdzdot-dxdzdotin*hnw-dzdzdotin*lm; 
dmdqm-dxndalsm*daldq-dxdqm*hm; 
dmdpui»dmdalsm*daldp-fdxdpm*hiT\; 

dmdthetom*dmdalBm*daldtheto-dxdthetom*hm+dzdtheto(n*lm; 
dmdalm-dmdalsm*dalda-dxdalm*hm ; 
dmdblm-dmdal sm* dal db - dxdblm* hm ; 

dndxdotm*rho*Ab* (ohm*R) A 2*R* (dcqsigdmu*dmudxdot+dcqsigdlamp*dlampdxdot) ; 
dndzdotm-rho*Ab* (ohm*R) A 2*R*dcqsigdlamp*dlampdzdot; 
dndrm--2/ohm*rho*Ab*(ohm*R) A 2*R*cqsig; 
dndthetom-rho*Ab*(ohm*R) A 2*R*dcqsigdtheto; 

% 

% return to CRUISE.M 
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0 . 


CTIt GRP.M 


% CTRGRP.M 
% CALLED BY CRUISE.M 

% Computes the stability derivatives of the tail rotor in cruise flight. 
% Computes the stability derivatives of the NOTAR in cruise flight. 

% Uses data loaded in the workspace by JANRAD.M and STAB.M. 

% 

if notar>*0 

% 

% FIX THIS 

% 

% Compute the stability derivatives of the tail rotor in cruise flight. 

% 

% dctsigdlamptxl.04 % table lookup, Prouty table 

9.5 

% dctsigdmut*-0.07 % table lookup, Prouty chap 3 

% dctsigdthetot*.659; % table 9.5 

dlampdydott=-l/((ohmt*Rt)*(l+dctsigdlampt*sigmat/(2*mut))); 
dydxdott=rho*Abt*(ohmt*Rt) A 2*dctsigdmut*dmudxdot; 
dydydott=rho*Abt* (ohmt*Rt) A 2*dctsigdlampt*dlampdydott ,- 
dydpt *dydydot t *ht; 
dydr t=-dydydot t *11; 

dydthetot*rho*Abt*(ohmt*Rt) A 2*dctsigdthetot; 

drdxdot t=dydxdot t * ht; 

drdydot t=dydydott *ht; 

drdpt»dydpt*ht ; 

drdrt=dydrt*ht; 

drdthetot=dydthetot*ht; 

dndxdott=-dydxdott*lt; 

dndydot t =-dydydot t *11; 

dndpt=-dydpt*lt ,- 

dndrt*-dydrt*lt; 

dndthetot=*dydthetot*lt; 

% 

% If tail rotor is to be used, zero out NOTAR derivatives: 

% 

dydxdotn=0;dydydotn=0;dydzdotn=0;dydpn=0dydrn=0;drdxdotn=0; 
drdydotn=0;drdzdotn=0;drdpn=0;drdrn=0;dndxdotn=0;dndydotn=0; 
dndzdotn=0; dndpn=0; dndm=0 ; 
dndphin=0;drdphin=0;dydphin=0; 
dnddelv=0;drddelv=0;dyddelv=0; 

elseif notar==l 

% 

% Compute the stability derivatives of the NOTAR in cruise flight. 

% 

qvin*.5*rho*vlmax A 2; 
dydxdotn=0; 

dydcsoar=qvin*bn A 2*max([(-V/125+1),0]); 

% downwash blown off circulation control section above 125 ft/sec (75 
kts) 

dcsoardcloar=l; 
dcloardpslot=2/pi; 

dpslotdcmu-dian/(2*bn)*sqrt(cmun*dian/(2*hslot))*1.5; 

dcmudvimax=-64*hslot/(dian*vlmax),- 

dvimaxdzdot*l; 

dcsoardswirl=-1/pi; 

dswirldydotal/vl; 
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dydydotnadydcsoar*dcsoardswirl*dswirldydot; 


dydzdotn>dydc8oar*dcaoardcloar*dcloardpslot*dpslotdanu*dcatudvimax*dvimax 
dzdot; 

dydpn»dydydotn*hn; 
dydm» - dydydotn*ln ; 
drdxdotn»dydxdotn*hn; 
drdydotn»dydydotn*hn; 
drdzdotn-dydzdotn*In; 
drdpn«dydpn*hn; 
drdrn*dydrn*hn ; 
dndxdotn--dydxdotn*In; 
dndydotn--dydydotn*In; 
dndzdotn*-dydzdotn*In; 
dndpns-dydpn*ln; 
dndm= - dydrn* In; 
dydphin«Ytmajcn/ (pi/2) ; 
drdphin*dydphin*htn; 
dndphin=-dydphin*1tn; 

% 

% If NOTAR is to be used, zero out tail rotor derivatives: 

% 

dydxdott=0;dydydott*0;dydpt*0;dydrt=0;dydthetot=0;drdxdott=0; 
drdydott=0;drdpt=0;drdrtaO;drdthetot=0;dndxdott=0;dndydott»0 ; 
dndpt * 0;dndrt * 0;dndthetot«0; 
mut=0;lampt=0;aot=0;blst=0;alst=0;locknot=0; 
drdthetot=0 ;dndthetot=0 ;dydthetot«=0 ; 
else 

disp(' ') 

disp(' ERROR IN CTRGRP.M') 

disp(' NO TAIL ROTOR OR NOTAR INSTALLED') 

disp(' CHECK INPUT DATA') 

disp{' ') 

disp(' Press any key to continue...') 

pause 

end 

% 

% return to CRUISE.M 
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H. CBODYGRP.M 

% CBODYGRP.M 
% CALLED BY CRUISE.M 

% Computes the stability derivatives of the fuselage, wing, verticle 
% fin and horizontal stabalizer in cruise flight. 

% Uses data loaded in the workspace by JANRAD.M STAB.M CRUISE.M 
% CMRGRP.M CTRGRP.M 

% 

» THIS SUBROUTINE MUST FOLLOW CMRGRP.M AND CTRGRP.M 

% 

% Compute the stability derivatives of che fuselage 

» 

dgamdzdotf * -1/V; 
dbetadydotf=l/V; 

detamfdxdotf*vfvl/(4*q*pi*R A 2)*(-dzdxdotm-2*T/V); 
detamfdzdot£=-vfvl/(4*q*pi*R A 2)*dzdzdotm; 
dalpfdxdotf=-detamfdxdotf; 
dalpfdzdotf«-(detamfdzdotf+dgamdzdotf); 

% 

% the folowing are from the curves in appendix A of Prouty: 

% and can be modified if necessary 

dfdalpf=-0.5; 

dlqdalpf=19.2; 

dsfqdbetaf=-55; 

dmqdalpf=445; 

dnqdbetaf=-205; 

drqdbetaf=58; 

% 

dxdxdotf=-2/v*f*q; 

dxdzdotf=-(-f*q-q*dfdalpf)*dalpfdzdotf; 
dydydotf=l/V*(q*dsfqdbetaf-f*q); 
dzdxdotf=0; 

dzdzdotf=(-f*q-q*dlqdalpf)*dalpfdzdotf; 
drdydot f=q*drqdbe t af *dbetadydotf; 
dmdxdotf=q*dmqdalpf*dalpfdxdotf; 
dmdzdotf=q*dmqdalpf*dalpfdzdotf; 
dndydotf=q*dnqdbetaf*dbetadydotf; 

% 

% Compute the stability derivatives of the wing 

% 

if wing==l 

dgamdzdotw=-l/V; 

detamhdxdotw=vwvl/(4*q*pi*R A 2)*{-dzdxdotm-2*T/V); 
detamhdzdotw=-vwvl/(4*q*pi*R A 2)*dzdzdotm;; 

detafhdzdotw=detafdalpfw*(1/(4*q*pi*R A 2)*dzdzdotm-dgamdzdotw); 
dalphdxdotw=-detamhdxdotw; 

dalphdzdotw=>(detamhdzdotw+detafhdzdotw+dgamdzdotw); 
dalphdzdbldotw=-detamhdzdotw*lw/V; 

\ 

dxdxdotw®-2 * cdow; 

dxdzdotw=q*Aw*aw* ((alpw-alplow) * (l-2*aw* (1+deliw) *Aw/pi/bw A 2) +alpw-iw) * 
dalphdzdotw; 

dxdzdbldotw=dxdzdotw*dalphdzdbldotw/dalphdzdotw; 

dzdxdotw=2/V*Zw-q*Aw*aw* (i+aw* (l+deliw) *Aw/pi/bw A 2* ( 2 * (alpw-alplow) * 
(alpw-iw)+(alpw-alplow) A 2)+cdow)*dalphdxdotw; 
dzdzdotw=-q*Aw*aw(l+aw*(l+deliw)*Aw/pi/bw A 2*(2*(alpw-alplow)*... 

(alpw-iw)+(alpw*alplow) A 2)+cdow)*dalphdzdotw; 
dzdzdbldotw=dzdzdotw*dalphdzdbldotw/dalphdzdotw; 
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% 

dzdqw*dzdzdotw*lw; 
dzdqw»dzdzdotw*lw; 

dmdxdotw*•dxdxdotw*hw+dzdxdotw*lw; 
dmdzdotw*•dxdzdotw*hw+dzdzdotw*lw; 
dxndzdbldotw* - dxdzdbldotw*hw+dzdzdbldotvr* lw ; 
dmdgw-dzdqw*lh; 

% 

drdrw«q*Aw*bw*bw*aw»alpw/4/V/V/Ixx; 

drdpws-q*Aw*bw*bw/4/V/V/lxx*pi*(l+3*ctw/crw)/(6*( 1 +ctw/crw)); 
else 
% 

% Zero stability derivatives of the wing when no wing is installed 

% 

dxebcdotwsO;dxdzdotw«0;dxdzdbldotw=0;dzdxdotw=0; 
dzdzdotwsO;dzdzdbldotwaO; 
dzdqw=0;dzdqw=0;dmdxdotw*0;dmdzdotw*o; 
dmdzdbldotw=0;dmdqw=0; 
drdrw=0;drdpw=0; 
end 
% 

% Compute the stability derivatives of the verticle fin 

% 

detafdbeta=.06; % assumed to be .06 because of little study of effect 
dbe tadydot v* 1 /V ; 

detatvdxdotv*-1/(4*qvq*q*Av)*(dydxdott-2*Tt/V); 
dalpvdxdotv=detatvdxdotv; 
detatvdydotv=dydydott/(4*qvq*q*At); 
detafdydotv=detafdbeta*dbetadydotv; 

dalpvdydotv=-(dbetadydotv+detatvdydotv+detafdydotv); 

% 

dxdxdotv=0; 
dxdydotv*0; 

dydxdotv=2/V*Yv+qvq*q*Av*av*dalpvdxdotv; 
dydydotv=l/(1-deldv/Yv)*... 

(qvq*q*A' ' ,av*dalpvdydotv+deldv* (-2/V+l/Tt*dydydott)) ; 

% 

dydpv=dydydotv*hv; 
dydrv=-dydydotv*1v; 
drdxdotv=dydxdotv*hv; 
drdydotv=dydydotv*hv; 
drdpv=dydpv*hv; 
drdrv=dydrv*hv ,• 
dndxdotv*-dydxdotv*1v; 
dndydotv=-dydydotv*1v; 
dndpv=-dydpv*1v; 
dndrv*-dydrv*lv; 

% 

dydde1v=q*qvq*Av* c1vertmax/delvmax; 
drddelv»dyddelv*hv; 
dnddelv=-dyddelv*1v; 

% 

% Compute the stability derivatives of the horizontal stabalizer 

% 

dgamdzdoth*-1/V; 

detamhdxdoth«vhvl/<4*q*pi*R A 2)*(-dzdxdotm-2*T*cos(altpp)/V); 
detamhdzdoth»-vhvl/(4*q*pi*R A 2)*dzdzdotm; 

detafhdzdothsdetafdalpfh*(1/(4*q*pi*R A 2)*dzdzdotm-dgamdzdoth); 
dalphdxdoth*-detamhdxdoth; 

dalphdzdoths-abs(detamhdzdoth+detafhdzdoth+dgamdzdoth); 







dalphdzdbldoth*-detamhdzdoth*lh/v; 

% 

etamh»vhvl*vl/V; %fix this 
eta£hsdeta£dalpfh*tho; 

Xh=Lhoriz+ein (tho- (etamh+etafh+gamc)) -Dhoriz*cos(tho- (etafh+etafh+gamc)) ; 
dxdxdoth=2/v*Xh+qhq*q*Ah*ah*((alph-alploh)*(l-2*ah*(l+delih)*Ah/pi/bh A 2) 
♦alph-ih)*... 
dalphdxdoth; 

dxdzdoth»qhq*q*Ah*ah*((alph-alploh)*(l-2*ah*(l+delih)*Ah/pi/bh A 2)+alph-i 
h)*... 

dalphdzdoth; 

dxdzdbidoth=cbcdzdoth*dalphdzdbldoth/dalphdzdoth; 

dzdxdoth*2/V*Zh-qhq*q*Ah*ah*(1+ah*(l+delih)*Ah/pi/bh A 2*(2+(alph-alploh)* 

(alph-ih)+(alph-alploh) A 2)+cdoh)*dalphdxdoth; 
dzdzdoth=-abs(-qhq*q*Ah*ah(1+ah* (l+delih) *Ah/pi/bh A 2* (2* (alph-alploh) *. . . 

(alph-ih)+(alph*alploh) A 2)+cdoh)*dalphdzdoth); 
dzdzdbldoth=dzdzdoth*dalphdzdbldoth/dalphdzdoth; 

% 

dzdqh=dzdzdoth+lh; 

dzdqh=dzdzdoth+lh; 

dmdxdoth=-dxdxdoth*hh+dzdxdoth*lh; 

dmdzdoths-dxdzdoth*hh+dzdzdoth*lh; 

dmdzdbldoth=-dxdzdbldoth*hh+dzdzdbldoth*lh; 

dmdqh=dzdqh*lh; 

% 

% return to CRUISE.M 
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Z. HOVER.M 

% HOVER.M 

% CALLED BY STABDER.M 

% Computes Che scability derivatives at a hover. 

% calls the FOLLOWING subroutines to compute stability derivatives 

% 

% HMRGRP 
% HTRGRP 
» TRIM 

« 

% computation of stability derivatives 

% the only derivatives important at hover are main and tail rotor 
format compact 

% evaluate dctsigdtheto dcqsigdtheto dctsigdthetot and dcqsigdthetot 

% 

vctsigszeros(1,2); 
vcqsigxzeros (1,2) 
vthetaoszeros(1,2); 
vctsigtxones(l,2); 
vcqsigt=ones(1,2); 
vthetaot=ones(1,2); 
eval(['load ',filenamel]) 
thetalx-twist; 

GW=GW*1.0005; % .5% higher 

disp(' ') 

disp(' *** FIRST OF TWO TRIM ITERATIONS ***') 

disp(' ') 
pause(3) 
trim 

% *** Calculation of output parameters from trim *** 

Qrotor=mean(DMpsi)*b; 
solidity=b*cblade/(pi*R); 

CQ=Qrotor/(Adisk*rho*Vtip A 2*R); 

CH*Hrotor/(Adisk*rho*Vtip A 2); 

CH_sig=CH/solidity; 

vctsig(l,2)=CT/solidity; 

vcqsig(l,2)xCQ/solidity; 

thetasave=thetao; 

vthe -=<- ' 1 , 2 ) xthetasave-.7*thetal; 

ohr;. .i ; 

sigm«i=j oiidity; 

g=32.174; 

A«pi*R*R; 
if notar==0 
Abt=bt*Rt*cot; 

At«pi*Rt*Rt; 
mut*Vinf/ohmt/Rt; 
sigmatxAbt/At ; 
lt=ltd-xcg; 
htxhtd-zcg; 
yt=ytd-ycg; 

locknotxrho*at*cot*Rt A 4/Ibt; 

Tt=Qrotor/lt; 

ctsigt-Tt/(rho*At*(ohmt*Rt) " 2 ); 
vlt=sqrt(Tt/2/rho/At); 

1ampt■-vlt/ohmt/Rt; 

aotx2/3*locknot*ctsigt/at-3/2*g*Rt A 2/(ohmt*Rt) A 2; 
thetaotx3/2*(4/at*ctsigt+sqrt(CQ/2))-3/4*thetalt; 

alStxO; 
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bl8t«0; 

theta75t-thetaot+.75*thetalt; 
vlt-sqrt(Tt/2/rho/At); 
aoat«theta75t+atan(vlt/ohmt/Rt/. 7) ; 
cdt».006+(aoat*at) A 2/pi/Rt/cot; 
phit«sqrt(ctsigt*sigmat/2); 
thet«4/at*ctsigt+sqrt(phit); 

cqsigt-sigmat/2*(at/2+(that-phit)*phit+cdt/4); 
and 

vctaigt(1,2)>ctsigt; 
vcqsigt(1,2).cqsigt; 
vthataot(1,2)-thetaot; 

% 

aave stabtemp vctsig vcqsig vthetao thetal filenamel vctaigt vcqsigt 

vthataot g 

clear 

load stabtemp 

aval(['load ',filenamel]) 

load stabtemp 

disp(' ') 

disp(' *** SECOND OF TWO TRIM ITERATIONS ***') 

disp(' ') 
pause(3) 
trim 

disp(' ') 

disp(' *** PERTURBATIONS COMPLETE ***') 

disp(' ') 

disp (' *** EVALUATING STABILITY DERIVATIVES ***') 

disp(' ') 

Qrotor=mean(DMpsi)*b; 
solidity=b*cblade/(pi*R); 

CQ=Qrotor/(Adisk*rho*Vtip A 2*R); 

CH=Hrotor/(Adisk*rho*Vtip A 2); 
vctsig(1,1)=CT/solidity; 
vcqsig(1,1)*CQ/solidity; 
thetasave vthetao; 
vthetao(1,1)*thetasave-.7*thetal; 

dctsigdtheto=(vctsig(1,1)-vctsig(1,2))/(vthetao(1,1)-vthetao(1,2)); 
dcqsigdtheto=(vcqsig(1,7'-vcqsig(1,2))/(vthetao(1,1)-vthetao(1,2)); 
ohmsomega ,- 
sigmassolidity; 

A=pi*R*R; 
if notar==0 

Abt=bt*Rt*cot; 

At»pi*Rt*Rt; 
mut=Vinf/ohmt/Rt; 
sigmat=Abt/At ; 
lt=ltd-xcg; 
ht*htd-zcg; 
yt*ytd-ycg; 

locknot«rho*at*cot*Rt A 4/Ibt; 

Tt*Qrotor/lt; 

ctsigt«Tt/(rho*At*(ohmt*Rt) A 2); 
vlt*sqrt(Tt/2/rho/At); 
lampt=-vlt/ohmt/Rt; 

aot=2/3*locknot*ctsigt/at-3/2*g*Rt A 2/(ohmt*Rt) A 2; 
thetaot*3/2*(4/at*ctsigt+sqrt(CQ/2))-3/4*thetalt; 
alst=0; 
bl8t*0; 

theta75t=thetaot+.75*thetalt; 
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vlt-sqrt(Tt/2/rho/At); 
aoat«theta75t+atan(vlt/ohmt/Rt/.7); 
cdt«.006+(aoat*at) A 2/pi/Rt/cot; 
phit-sqrt(ctsigt*sigmat/2); 
thet*4/at*ctsigt+sqrt(phit); 

cqsigt*sigmat/2*(at/2*(thet-phit)*phit+cdt/4); 
end 

vctsigt(1,1)sctsigt; 
vcqsigt(1,1)«cqsigt; 
vthetaot(1,1)=thetaot; 

dctsigdthetots(vctsigt(1,2)-vctsigt(1,1))/(vthetaot(1,2)-vthetaot(1,1)); 
dcqsigdthetot-(vcqsigt(1,2)-vcqsigt(1,1))/(vthetaot(1,2)-vthetaot(1,1)); 
%%%%%%%%»%%%%%%%%%%%%%%%t%%%%%t%%%%%%»%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% CONFIGURATION CALCULATIONS 
% converting data from crank 
theta7«thetao; 
thetao«theta7-.7*thetal; 
ctsig«CT/solidity; 
chsig*CH/solidity; 
cqsig*CQ/solidity; 
altppc-alphaT; 

Al*-thetalc; 

Bl«-thetals; 
ao*betao; 
vl«mean(vi); 

V.Vinf; 

% Main rotor 

9-32.2; 

hm*hmd-zcg; 

ym=ymd-ycg; 

lnudmd-xcg; 

c-cblade; 

Ab=c*R*4; 
lamp=0; 
mu»0 ; 

Ct=CT; 
f =Af v; 

thetat=4/a*ctsig+sqrt(sigma*ctsig/2); 
lockno=rho*a*c*R A 4/lb; 
theta75*thetao+.75*thetal; 
numl*thetao* (8/3+32/45*mu A 3/pi) ,- 
num2*thetal*(2+mu A 4/12); 
num3«lamp* (2-mu*mu/2) ; 
denl*l+3/2*mu*mu-5*inu A 4/24 ; 
als=mu* (numl+num2+num3) /denl-Bl ; 
bls-als*12/lockno*e/R/(l+e/R/3); 
m*GW/g; 

Ic»Ixx*Izz; 
if notar--l; 

% NOTAR 

vlmax«1.2*vl; % downwash at the boom slots (NOTAR) 

delvmax>cpi/4 ; 

ddelvddelp=delvmax/maxr; 

ltn-lttnd-xcg; 

lnsltnd-xcg; 

yn*ytnd-xcg; 

hn«htnd-zcg; 

htn*hn; 

Tt*cqsig*sigma*rho*A*(ohm*R) A 2*R/ltn; 
cmun*.55; % typical value 
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bn*.5*R; 

hslot*.028*dian/2; 
am-bn/dian; 
end 

% trim conditions 
tho*-T*(im+als)/GW; 
pho--(T*blS+Tt)/GW; 
wo*0; 
vo*0; 

UO*0; 

%%%%%%%»%%%%%%%%%%%%%%%%%%%%%%%%%«%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% STABILITY CALCULATION 

% 

hmrgrp 

htrgrp 

% 

% computation of A,B,C,D matrices 

% 

All* (dxdxdotm)/m; 

A12* (dxdxdotm)/m; 

A13* (dxdqm)/m-wo; 

A14* -g*cos(tho); 

A15* (dxdydotm)/m; 

A16* (dxdpm) /m; 

A17* 0; 

A18* (dxdrm)/m+vo; 

% 

A21* (dxdxdotm)/m; 

A22* (dxdxdotm)/m; 

A23* (dxdqm)/m+uo; 

A24* -g*cos(pho)*sin(tho); 

A25* (dxdydotm+dxdydott)/m; 

A26* (dxd^pm) /m-vo; 

A27* -g*sin(pho)*cos(tho); 

A28* (dxdrm)/m; 

% 

A31* (dmdxdotm)/Iyy; 

A32* (dmdxdotm)/Iyy; 

A33* (dmdqm)/Iyy; 

A34* 0; 

A3S* (dmdydotm+dmdydott+dmdydotn)/Iyy; 

A36* (dmdpm)/Iyy; 

A37* 0; 

A38= (dmdrm+dmdrt+dmdm) /Iyy; 

% 

A41* 0; 

A42* 0; 

A43* cos(pho); 

A44* 0; 

A45- 0; 

A46* 0; 

A47* 0; 

A48* -sin(pho); 

% 

A51* (dydxdotm)/m; 

A52* (dydxdotm+dydzdotn) /m; 

A53* (dydqm)/m; 

A54* -g*sin(pho)*sin(tho); 

A55* (dydydotm+dydydott+dydydotn)/m; 

A56* (dydpm+dydpt+dyc^n)/m+wo; 
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A57* g*cos(pho)*cos(tho); 

A58- (dydrm+dydrt+dydm) /m-uo; 

% 

A61« (Izz* (drdxdotm+drdxdotn) +Ixz* (dndxdotm+dndxdotn) ) /Ic; 

A62» (Izz* (drdzdotm+drdzdotn) +Ix»* (dndxdotm+dndxdotn)) /Ic; 

A63» (Izz* (drdtjn) +Ixz* (dndqm) )/Ic; 

A64- 0; 

A 6 5 

(Izz*(drdydotm+drdydot t+drdydotn)+Ixz*(dndydotm+dndydot t+drdydotn))/Ic; 
A66» (Izz*(drdpm+drdpt+drdpn)+lxz*(dndpm+dndpt+dndpn))/Ic; 

A67. 0; 

A68> (Izz* (drdrm+drdrt+drdrn) +Ixz* (dndnn+dndrt+dndm) ) /Ic; 

% 


A71- 

A72«= 

A73« 

A74» 

A75- 

A76- 

A77« 

A78> 

t 

A81- 

A82« 

A83- 

A84* 

A 


0 
0 

sin(pho)*tan(tho) 
0 ; 

0 
1 
0 

cos (pho) *taui (tho) 


(Ixz*(drdxdotm+drdxdotn)+Ixx*(dndxdotm+dndxdotn))/Ic; 
(Ixz*(drdzdotm+drdzdotn)+Ixx*(dndzdotm+dndzdotn))/Ic; 
(Ixz*(drdqm)+ ixx* (dndqm))/Ic; 

0 ; 

8 5 


(Ixz*(drdydotm+drdydott+drdydotn)+Ixx*(dndydotm+dndydott+dndydotn))/Ic; 
A86* (Ixz* (drd^m-«-drdpt>dr<^}n) +Ixx* (dnc^im+dndpt+dndpn) ) /Ic; 

A87» 0; 

A88* (Ixz*(drdrm+drdrt+drdrn)+Ixx*(dndrm+dndrt+dndm))/Ic; 

% 

% longitudinal plant augmented X*(u w q theta] 

Flonaug*> [All A12 A13 A14;A21 A22 A23 A24/A31 A32 A33 A34 ; 

A41 A42 A43 A44] ; 

Plonaug-poly(Flonaug); 

Rlonaug«roots(Plonaug); 


% 


% Lateral plant augmented with X-[v p phi r psi]' 
Flataugs [AS5 A56 A57 A58 0;A65 A66 A67 A68 0; 

A75 A76 A77 A78 0;A85 A86 A87 A88 0;0 0 0 1 0] ; 
Plataug=poly(Flataug); 

Rlataug«roots(Plataug); 

% coupled plant 

Amat« [All A12 A13 A14 A15 A16 A17 A18 0; 

A21 A22 A23 A24 A25 A26 A27 A28 0; 

A31 A32 A33 A34 A35 A36 A37 A38 0; 

A41 A42 A43 A44 A45 A46 A47 A48 0; 

A51 AS2 A53 A54 ASS AS6 A57 A58 0; 

A61 A62 A63 A64 A6S A66 A67 A68 0; 

A71 A72 A73 A74 A75 A76 A77 A78 0; 

A81 A82 A83 A84 A85 A86 A87 A88 0; 

00000001 0 ]; 


Pcoup*poly (Amat) ; 
Rcoup»roots(Pcoup); 

% 


Bll* dxdblm*dblmddele/m; 

B12« dxdthetom*dthetomddelc/m; 
B13« dxdalm*dalmddela/m; 

B14« 0; 
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% 

B21» dzdblm*dblmddele/m; 

B22- dzdthetom*dthetomddelc/m; 

B23- 0; 

B24- 0; 

% 

B31- dmdblm*dblmddele/Iyy; 

B32- dmdthetom*dthetomddelc/Iyy; 

B33» dmdalm*dalmddela/Iyy; 

B34« (dmdphin*dphinddelp+dmdthetot*dthetotddelp)/lyy; 

% 

B41- 0; 

B42- 0; 

B43- 0; 

B44- 0; 

% 

B51« dydblm*dblmdde1 e /m ; 

B52. dydthetom*dthetomddelc/m; 

B53» dydalm*dalmddela/m; 

B54« (dydphin*cphinddelp+dydthetot*dthetotddelp) /m; 

% 

B61» Izz*drdblm*dblmddele/Ic; 

B62 = (Izz*drdthetom*dthetomddelc+Ixz*dndthetom*dtheto«nddelc)/Ic; 
B63= Izz*drdalm*dalmddela/Ic; 

B64» (Izz* (drdphin*dphinddelp+drdthetot*dthetotddelp) +. . . 

Ixz*(dndphin*dphinddelp+dndthetot*dthetotddelp))/Ic; 

* 

B71* 0; 

B72* 0; 

B73* 0; 

B74- 0; 

% 

B81= Ixz*drdblm*dblmddele/Ic; 

B82= (Ixz*drdthetom*dthetomddelc+Ixx*dndthetom*dthetomddelc)/Ic; 
B83® Ixz*drdalm*dalmddela/Ic; 

B84= (Ixz*(drdphin*dphinddelp+drdthetot*dthetotddelp)+... 

Ixx # (dndphin*dphinddelp+dndthetot*dthetotddelp))/Ic; 

% 

Glonaug-[Bll B12 B13 B14; 

B21 B22 B23 B24; 

B31 B32 B33 B34; 

B41 B42 B43 B44]; 

% 

Glataug*[B51 BS2 B53 B54; 

B61 B62 B63 B64; 

B71 B72 B73 B74; 

B81 B82 B83 B84; 

0 0 0 0 ]; 

% coupled input matrix 

Bmat«[Bll B12 B13 B14;B21 B22 B23 B24; 

B31 B32 B33 B34.-B41 B42 B43 B44; 

B51 B52 B53 B54.B61 B62 B63 B64 ; 

B71 B72 B73 B74;B81 B82 B83 B84;0 000]; 

» 

% cross coupling 
xcouplemi2/lockno*e/R/(l+e/3/R); 

% designed damping 
desdmdq«dmdqm; 
desdrdp-drdpm+drdpt+drdpn; 
desdndr*dndrm+dndrt+dndrn; 












% now cooper harper pilot rating 
prpitch-desdmdq/Iyy; 
prroll* (drdpm+drdpt+drdpn) /Ixx; 
pryaw*desdndr/Izz; 

% control power 
cppitch*B31*lyy; 
cprol1-B63*Ixx; 
cpyaw>B84*Izz; 
cpipitch«B31; 
cpiroll«B63; 

Cpiyaw»B84; 

% 

thetao«theta7; 
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J. 


HMRQRP.M 


t HMRGRP.M 
t CALLED BY HOVER.M 

% Computes the basic main rotor derivatives at a hover 
% Computes the stability derivatives of the main rotor at a hover 
% Uses data loaded in the workspace by JANRAD.M and STAB.M 

% 

% Compute the basic mainrotor derivatives at a hover 

% 

dmudxdot»1/ohm/R; 
dlanpdzdotsdmudxdot; 

daldmu»8/3*thetao+2*thetal-2*vl/ohm/R; 

dbldmu»4/3*ao; 

dctsigdlamp*inv(8/a+(sqrt(sigma/2)/(sqrt(ctsig)))); 
dcqsigdlamps-a/4*(theta75-2*vl/ohm/R); 

daldq=-16/(lockno*ohm*(1-e/R)*2)-12*e/R/(lockno*ohm*(1-e/R)*3); 
dbldp«daldq; 

daldp»l/ohm* (1- (192*e/R/ (lockno A 2* (1-e/R) '“'5))) ; 
dbldq=-daldp; 

dalda=12*e/R/(lockno*(1-e/R)“3); 
dbldb=dalda; 

daldb=-l/(1+((144*(e/R) A 2)/(lockno A 2*(1-e/R)‘6))); 
dblda=-daldb; 

dchsigda*3/2*ctsig*(l-a/18*theta75/ctsig); 
dcysigdb=dchsigda; 

dmdals=3/4*e/R*Ab*rho*R*(ohm*R) A 2*a/lockno; 
drdbls=dmdals; 

% 

% Compute the mainrotor stability derivatives at a hover 

% 

dxdxdotm=-rho*Ab*(ohm*R) A 2*dchsigda*daldmu*dmudxdot; 
dxdydotm=-rho*Ab*(ohm*R) A 2*dcysigdb*dbldmu*dmudxdot; 
dxdzdotm= 0; 

dxdqm=-rho*Ab*(ohm*R) A 2*dchsigda*daldq; 
dxdpm*-rho*Ab*(ohm*R) A 2*dchsigda*daldp ; 
dxdrm=0; 

dxdthetom=-rho*Ab*(ohm*R) A 2*(als+im)*dctsigdtheto; 

dxdalm=-rho*Ab*(ohm*R) A 2*dchsigda*dalda; 

dxdblm=-rho*Ab*(ohm*R) A 2*dchsigda*daldb; 

dydxdotm=rho*Ab*(ohm*R) A 2*dcysigdb*dbldmu*dmudxdot; 

dydydotm=dxdxdotm; 

dydzdotm=0; 

dydqm=dxdpm; 

dydrm=0; 

dydpm=-dxdqm; 

dydthetom«rho*Ab*(ohm*R) A 2*(2*bls)*dctsigdtheto; 

dydalm=dxdblm; 

dydblm=-dxdalm; 

dzdxdotm-0; 

dzdydotm=0; 

dzdzdotm=-rho*Ab*(ohm*R) A 2*dctsigdlamp*dlampdzdot; 
dzdpm=0; 
dzdqm=0; 
dzdrm=0; 

dzdthetom=-rho*Ab*(ohm*R) A 2*dctsigdtheto; 
dzdblm=0; 

drdxdotm«drdbls*dbldmu*dmudxdot+dydxdotm*hm+dzdxdotm*hm; 
drdydotm*-dmdals*daldmu*dmudxdot+dxdxdotm*hm; 
drdzdotm=dzdzdotm*ym; 








drdqm*drdblB*dbldq-t-dydqm*hni ; 
drdrm*0; 

drdpm-drdblB*dbldp-*-dydpm*hin; 

drdthetom»dydth«tom*hm+dzdtheto«n*ym ; 

drdalm»drdbl8*dblda+dydalra*hm; 

drdblm-drdbl•* dbldb+dydblm*hn» ; 

dmdxdotm*drodal»*daldBSU*dmudxdot-dxdxdotm*hin; 

dmdydotm-drdxdotm ; 

dmdzdotm-dz dzdotm*lm; 

dmdqm-dmdal s *daldq• dxdqpi*hm 

dmdrmaO; 

dmdpm*dmdals*daldp-dxdpm»hin; 
dmdthe torn* - dxdtheto*n*hm+dzdthetom* lm ; 
dmdalm«dmdal8*dalda-dxdalm*hm; 
dmdblm«dmdalB*daldb-dxdblm*hm; 
dndxdotm«0; 

dndzdotm«rho*Ab* (ohm*R) ' ,l 2*R*dcqsigdlainp*dlainpdzdot ; 
dndydotn.-O ; 
dndqgmaO; 
dndpm^O; 

dndnn«-2*rho*Ab*(ohm*R) A 2*R*cqsig; 
dndthetom»rho*Ab*(ohm*R) A 2*R*dcqsigdtheto; 

% 

% return to HOVER.M 
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X. 


HTRGRP. M 


% HTRGRP.M 
% CALLED BY HOVER.M 

% Computes the basic tail rotor or NOTAR derivatives at a hover 
% Computes the stability derivatives of the tail rotor or NOTAR at a hover 
t OR 

% Computes the stability derivatives of the NOTAR at a hover 

% 

% Uses data loaded in the workspace by JANRAD.M and STAB.M 

* 

if notar*sO 

% 

% Compute the basic tail rotor derivatives at a hover. 

% 

dmudxdott«l/ohmt/Rt; 
dlampdydotts-dmudxdott; 

daldmut=8/3*thetao+2*thetal-2*vl/ohmt/Rt; 

dbldmut*4/3*aot; 

dctsigdlampt=inv(8/at+(sqrt(sigmat/2)/(sqrt(ctsig)))); 

dcqsigdlampt=-at/4*(theta75t-2*vlt/ohmt/Rt); 

daldqt«0; 

dbldpt=0; 

daldptsO; 

dbldqtsO; 

daldat=0; 

dbldbtsO; 

daldbtsO; 

dbldatsO; 

dchsigdat=3/2*ctsig*(1-at/18*theta75t/ctsig); 
dcysigdbt=dchsigdat; 
dmdalst=0; 
drdblst=0; 

% 

% Compute the tail rotor stabitity derivatives at a hover 

% 

dydydott=rho*Abt*(ohmt*Rt)“2*dctsigdlampt*dlampdydott; 
dydpt=dydydott*ht; 
dydrt=-dydydott *11; 

dydthetot=rho*Abt* (ohmt*Rt) '*2*dctsigdthetot ; 

drdydott=dydydott*ht; 

drdrt=dydrt*ht; 

drdpt=dydpt*ht; 

drdthetot=dydthetot*ht; 

dmdydott=rho*Abt* (ohmt*Rt) '‘2*Rt*dcqsigdlampt*dlampdydott; 
dmdrt =-dmdydott *11; 

dmdthetot=rho*Abt*(ohmt*Rt)*2*Rt*dcqsigdthetot; 

dndydott*-dydydott *11; 

dndpt=-dydpt*lt; 

dndrt=-dydrt*lt; 

dndthetot=-dydthetot*lt; 

% 

% If tail rotor is used, zero out NOTAR derivatives: 

% 

dydxdotn=0; dydydotn=0; dydzdotn=0; dydpn=0; dydm=0 ; 
dzdydottsO; 

drdxdotnsO; drdydotn=0; drdzdotn=0; drd^n=0; drdm=0 ; 
dmdydotnsO; dmdm=0 ; 

dndxdotn=0;dndydotn=0;dndzdotn=0;dndpn=0;dndrn*0; 
dmdphin=0;drdphin=0;dndphin»0;dydphin=0; 
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elseif notar««l 

% Compute the NOTAR derivatives at a hover. 

% 

qvin-.5*rho*vlmax A 2; 
dydxdotrt»0; 
dydcsoar«qvin*bn A 2; 
dcsoardcloar-1; 
dcloardpslot»2/pi ; 

dpslotdcmu>dian/(2*bn)*sqrt(cmun*dian/(2*hslot))*1.5; 

dcmudvimax--64*hslot/(dian*vlmax); 

dvimaxdzdotsl; 

dcsoardswir1■-1/pi; 

dswirldydot*l/vi; 

dydydotn»dydcsoar*dcsoard8wirl*dswirldydot; 
dydzdotn*dydcsoar*dcsoardcloar*dcloardpElot*dpslotdcmu*dcmudviinax*dviniax 

dzdot; 

dydpn*dydydotn*hn; 
dydm- - dydydotn* In ; 
drdxdotn«dydxdotn*hn; 
drdydotn»dydydotn * hn; 
drdzdotn-dydzdotn*In; 
drdpn=dydpn*hn; 
drdm«dydm*hn ; 
dmdydotn*sO ; 
dmdm= - dmdydotn* In ; 
dndxdotn--dydxdotn*In; 
dndydotns-dydydotn*In; 
dndzdotn=-dydzdotn*ln; 
dndpn=-dydpn*In; 
dndm= - dydra* In ; 
dydphin-Ytmaxn/(pi/2); 
drdphin*dydphin*htn; 
dmdphin=0; 

dndphin*-dydphin*1tn; 

% 

% If NOTAR is to be used, zero out tail rotor derivatives: 

% 

dydydott=0dydpt=0; dydrt=0; dydthetot=0; drdydott=0 ; 
dzdydott«0; 

drdrt=0;drdpt«0;drdthetot=0; 
dmdydottsO;dmdrt*0;dmdthetot»0; 
dndydott*0dndpt=0 ;dndrt=0 ;dndthetot*0; 
lampt-0;mut=0;aot=0;blst=0;alst=0;locknot=0; 
dmdthetot«0;dndthetot=0;drdthetot=0;dydthetot*0; 
else 

disp(' ') 

disp(' ERROR IN HTRGRP.M') 

disp(' NO TAIL ROTOR OR NOTAR INSTALLED') 

disp(' CHECK INPUT DATA') 

disp(' ') 

disp(' Press any key to continue...') 
pause 
end 
% 

% return to HOVER.M 
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L. 


STABOOT.H 


% STABOUT.M 

% 

% calls 
% CMDBWPLH 
% CMDBWPLC 

% 

♦ *** Stability and Control Output Subroutine *** 

pack 

format compact 
clc 

disp{' ') 
disp(' ') 

disp('Do you want the results displayed on screen?') 
disp(' ') 

disp(' NOTE: if you want a hard copy of the plots, you must') 
disp(' select (l) and view them on the screen first.') 
flag=l; 

answer=input('1. yes 2. no >> '); 
while flag>0 

if answer == 2, 
flag=0; 

elseif answer == l, 

% 

% output to screen 

% 

clc 

disp(' ') 

disp(' *** STABILITY AND CONTROL PROGRAM ***') 

disp(' *** SCREEN VIEW MENU ***') 

disp(' ') 

disp('What do you want to see?') 
dispC ') 

disp('1. Input data.') 
disp('2. Calculated data.') 
disp(‘3. State Matrices.') 

disp('4. Eigenvalues of the plants and plots of the roots.') 
disp('5. Key control parameters.') 
disp('6. Open loop transfer plots.') 
disp(' ') 

disp('0. Exit screen view.') 
disp(' ') 

choice*input('Enter a number; '); 
if choice**!. 


clc 

disp(' ') 
disp(' 

eval (['disp('' 
disp(' ') 
disp (' 
disp(' ') 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
disp(' ') 
disp(' 
disp(' ') 
fprintf(' 


*** INPUT DATA (screen 1 of 8) ***') 

',filenamel,''')']) 

Flight Conditions') 

Forward velocity * %6.0f kts\n',Vinf/I.69) 
Temperature * %6.0f degs F\n',tenp) 
Pressure altitude * %6.0f ft\n',PA) 
Auxiliary thrust * %6.0f lbs\n',Taux) 

Fuselage') 

Gross weight = %6.0f lbs\n',GW) 
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fprintf(' 
fprintf(* 
fprintf(' 
fprintf(' 
fprintf('CG 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
disp(' ') 
disp('press 
pause 
clc 

disp(' ') 
disp(' 
eval(['disp 
disp(' ') 
disp (’ 
disp(' ') 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
if airfoil== 
disp(' 
else 

disp(' 


Equivalent flat plate area 
Vertical projected area 
CG height above waterline 
CG fuselage station 
position right of buttline 

Ixx 

lyy 

Izz 

Xxz 

Downwash ratio 
any key to continue. . .') 


*6.If ft A 2\n',Afh) 

%6.If ft A 2\n',Afv) 

%6.If ft\n',zcg) 

%6.If ft\n',xcg) 

%6. If ft\n',ycg) 

%6. if slug ft A 2\n',lxx) 
%6.1f slug ft A 2\n',Iyy) 
%6.If slug ft A 2\n',Izz) 
%6.1f slug ft A 2\n',Ixz) 
%6.2f \n',vfvl) 


*** INPUT DATA CONTINUED (screen 2 of 8) ***’) 
('' ',filenamel, ) 

Main Rotor') 


1 , 


Number of blades 
Rotor radius 
Blade chord 
Blade twist 


%6.Of \n',b) 

%6.If ft\n',R) 

%6.1f ft\n',cblade) 

%6.2f degs\n',twist*57.3) 


Blade airfoil 


HH-02') 


Blade airfoil 


VR-12') 


end 

fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
disp(' ') 
disp(' ') 
disp('press 
pause 
clc 


Blade lift curve slope 
Blade weight 
Rotational velocity 
Blade grip length 
Hinge offset 
Flapping moment of inertia 
Hub height above waterline 
Hub fuselage station 
Hub position rt of buttline 
Mast incidence 


any key to continue...') 


%6.2f \n',a) 

%6.If lbs\n',wblade) 
%6.2f rads/sec\n',omega) 
%6.If ft\n',grip) 

%6.If ft\n',e) 

%6.If slug ft A 2\n',Ib) 

%6.If ft\n',hmd) 

%6.If ft\n',lmd) 

%6.1f ft\n',ymd) 

%6.2f deg\n',im*57.3) 


disp(' ') 
disp(' 

eval (['disp( 
disp(' ') 
disp(' 
disp(' ') 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 


*** INPUT DATA CONTINUED (screen 3 of 8) ***') 

',filenamel,'*')']) 

Tail rotor (zeros if using NOTAR)') 


Number of blades = 
Blade chord = 
Blade radius = 
Lift curve slope = 
Rotational velocity = 
Flapping moment of inertia » 
Delta-3 angle s 
Blade twist « 
Hub height above waterline s 
Hub fuselage station e 


*6.If \n',bt) 

%6.If ft\n',cot) 

%6.If ft\n',Rt) 

*6.2f \n',at) 

V6.2f rad/sec\n*,ohmt) 

%6.If slug ft A 2\n',Ibt) 
%6.2f deg\n',delta3*57.3) 
%6.2f deg\n',thetalt*57.3) 
%6. If ft\n',htd) 

*6.If ft\n',ltd) 
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fprintf(' Hub position rt of buttline * %6.1f ft\n',ytd) 
disp(' ’) 

disp('press any key to continue...') 

pause 

clc 

disp(' ') 

disp(' *** INPUT DATA CONTINUED (screen 4 of 8) * 

eval (['disp('' filenamel,) 

dispj' ') 

disp(' NOTAR (zeros if using tail rotor)') 


%6 . If ft A 2\n',htnd) 

¥6.If ft A 2\n',ltnd) 

%6 . If ft A 2\n',ytnd) 

%6 .If ft A 2\n',dian) 

%6.2f deg\n',swirl*57.3) 
%6.If lbs\n',Ytmaxn) 


disp(' NOTAR (zeros if using tail rotor)') 

disp(' ') 

fprintf(' Height above waterline * %6.1f ft A 2\n', 

fprintf(' Fuselage station * ¥6.If ft A 2\n', 

fprintf!' Position right of buttline » %6.1f ft A 2\n', 
fprintf!' NOTAR boom diameter * ¥6.If ft A 2\n', 

fprintf(' Swirl angle at boom * ¥6.2f deg\n',£ 

fprintf(' Maximum thruster force = ¥6. If lbs\n',)i 

fprintf!' Thrust fuselage station = %6.1f ft A 2\n', 

disp(' ') 

disp('press any key to continue...') 

pause 

clc 

disp (' ') 

disp(' *** INPUT DATA CONTINUED (screen 5 of 8) 

eval! ['disp('' '.filenamel,*'')' 

disp(' ') 

disp(' Wing') 


Thrust fuselage station = %6.1f ft A 2\n',lttnd) 


disp(' 
disp(' ') 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
disp(' ') 


Area 
Span 
CL 
CDo 
Tip cord 
Root cord 
wing efficiency factor 
Zero lift angle 
Angle of incidence 
Lift curve slope 
Height above waterline 
Fuselage station 
Position right of buttline 
Rotor downwash ratio 
Fuselage downwash ratio 


¥6.If ft A 2\n',Swing) 

¥6.If ft\n',bwing) 

¥6.2f \n',CLwing) 

¥6.4f \n',CDowing) 

¥6.If ft\n',ctw) 

¥6.If ft\n',crw) 

¥6.2f \n',ewing) 

¥6.2f deg\n',alplow*57.3) 
¥6.2f deg\n',iw*57.3) 
¥6.2f \n',aw) 

¥6.If £t\n',hwd) 

¥6.If ft\n',lwd) 

¥6.If ft\n',ywd) 

¥6.2f \n',vwvl) 

¥6.2f \n',detafdalpfw) 


disp(' ') 

disp('press any key to continue...') 
pause 
clc 

disp(' ') 

disp(' *** INPUT DATA CONTINUED (screen 6 of 8) * 

eval(['disp('' ',filenamel,'*')')) 

disp(' ') 

disp(' Horizontal tail') 

disp(' ') 

fprintf(' Area * ¥6.If ft A 2\n',Shoriz) 

fprintf(' Span * ¥6.If ft\n',bhoriz) 

fprintf(' CL = ¥6.2f \n',CLhoriz) 

fprintf(' CDo ■ ¥6.4f \n',CDohoriz) 

fprintf(' Zero lift angle * ¥6.2f deg\n',alploh*57.3) 

fprintf(' Angle of incidence * ¥6.2f deg\n',iw*57.3) 

fprintf(’ Lift curve slope = ¥6.If \n',ah) 




fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
disp(' ') 
disp('press 
pause 
clc 

disp(' ') 
disp(' 

eval ( ['disp(' 
disp(' ') 
disp(' 
disp(' ') 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf <' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
disp(' ') 
disp('press 
pause 
clc 

disp(' ') 
disp(' 

eval(['disp(' 
disp(' ') 
disp(' 
disp(' ') 
fprintf(' 
fprintf(' 
fprintf( 


Height above waterline 
Fuselage station 
Position right of buttline 
Dynamic pressure ratio 
Rotor downwash ratio 
Fuselage downwash ratio 


¥6.If ft\n',hhd) 

%6. If ft\n\lhd) 

¥6.If ft\n',yhd) 

¥6.2f \n',qhq) 

%6.2f \n',vhvl) 

%6.2f \n',detafdalpfh) 


any key to continue...') 


INPUT DATA CONTINUED (screen 7 of 8) *' 

*,filenamel,**')']) 


Vertical tail') 


Area 

Span 

CL 

CDo 

Height above waterline 
Fuselage station 
Position right of buttline 
Zero lift angle 
Maximum Cl 
Dynamic pressure ratio 
Lift curve slope 


%€. If 
%6. If 
%6.2f 
¥6.4f 
¥6. If 
¥6. If 
¥6. If 
¥6.2f 
¥6. If 


¥ 6 . 

¥ 6 . 


2f 

2f 


ft A 2\n',Svert) 
ft\n',bvert) 

\n',CLvert) 

\n',CDovert) 

ft\n*,hvd) 

ft\n',lvd) 

ft\n',yvd) 

deg\n',alplov*57.3) 

\n',clvertmax) 

\n', qvq) 

\n',av) 


any key to continue. 


) 


INPUT DATA CONTINUED (screen 8 of 

',filenamel,''' 

Rigging') 


8 ) * 
)']) 


Long 

Lat 


cyclic pitch/inch defl = ¥6.2f deg/in\n',dblmddele*57.3) 
cyclic pitch/inch defl = ¥6.2f deg/in\n',dalmddela*57.3) 
Collective pitch/inch defl = ¥6.2f 

deg/in\n',dthetomddelc*57.3) 
fprintf(' Tail rotor 
deg/unit\n',dthetotddelp*57.3) 


p. :ch 


change/defl 


¥ 6 . 2 f 


Max deflection of control') 
from neutral for NOTAR = 
Displacement of anti-torque') 
control until full rudder = 


disp(' 
fprintf(' 
disp (' 
fprintf(' 
disp(' ') 

disp('press any key to continue...') 
pause 

elseif choice«2, ¥ calculated data 


¥6.2f units\n',dphinddelp*57, 
¥6.2f units\n',maxr) 


3) 


clc 

disp(' ') 
disp(' 

eval (['disp('' 
disp(' ') 
disp(' 
disp(' ') 
fprintf(' 
fprintf(' 


*** CALCULATED DATA (screen 1 of 2) ***') 

',filenamel,''')']) 

Main Rotor') 

thrust = ¥6.If lbs\n',T) 
torque = ¥6.If ft-lbs\n',Qrotor) 
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advance ratio 
inflow parameter wrt TPP 
Tip path angle 
Rotor coning angle 
1st lat cyclic term-Al 
1st long cyclic term-Bl 
lateral flapping 
longitudinal flapping 
Lock number 


’) 


CALCULATED 


fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf{' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
dispC ') 

disp('press any key to continue., 
pause 
clc 

disp(' ') 
disp(' 

eval ( ['disp(' 
disp(' ') 
disp(' 
disp(' ') 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
dispr ') 
disp(' ') 

disp('press any key to continue 
pause 

elseif choice«*3, % state matrices 


%6 

%6 

%6 

%6 

%6 

%6 

%6 

%6 

V6 


. If \n*,mu) 

. 3f \n* , lamp) 

.If degs\n',altpp*S7.3) 
.If degs\n',ao*57.3) 

.If degs\n',A1*57.3) 

.If degs\n',B1*57.3) 

.2f degs\n',bls*57.3) 

.2f degs\n',als*57.3) 
.If \n',lockno) 


DATA (screen 
',filenamel, 


2 of 2)***') 
'")']) 


Tail Rotor (zero if NOTAR)') 


tail rotor thrust = %6 
advance ratio • %6 
inflow parameter * V6 
Rotor coning angle ■ %6 
lateral flapping « %6 
longitudinal flapping * %6 
Lock number « V6 


) 


.If lbs\n',Tt) 

. If \n',mut) 

.3f \n',lampt) 

.If degs\n',aot*57.3) 

.2f degs\n'.blst*^.3) 
.2f degs\n',ai8t*57.3) 
.If \n',locknot) 


disp(' ') 

disp('Longitudinal uncoupled plant (A or F depending on notation)') 

disp('States are [u w q theta]') 

disp(' ') 

disp(Flonaug) 

disp(' ') 

disp(' ') 

disp ('Longitudinal uncoupled input matrix (B or G depending on notation)') 

disp('Inputs are [longitudinal cyclic, collective, lateral cyclic, 

pedals]') 

disp{' ') 

disp(Glonaug) 

disp(' ') 

disp('press any key to continue...') 

pause 

clc 

disp(' ') 

disp ('Lateral/directional uncoupled plant (A or F depending on notation)') 

disp('States are [v p phi r psi]') 

disp(' ') 

disp(Flataug) 

disp(' ') 

disp (' ') 

disp('Lateral/directional uncoupled input matrix (B or G depending') 
disp('on notation)') 

disp('Inputs are [longitudinal cyclic, collective, lateral cyclic, 
pedals]') 
disp(Glataug) 
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disp(' ') 

disp('press any key to continue...') 

pause 

clc 

disp(' ') 

disp{'Coupled plant (A or F depending on notation)') 
disp{'States are [u w q theta v p phi r psi]') 
disp(' ') 
disp(Amat) 

disp('press any key to continue...') 

pause 

clc 

disp(' ') 

disp('Coupled input matrix (B or G depending on notation)') 

disp('Inputs are [longitudinal cyclic, collective, lateral cyclic, 

pedals]') 

disp(' ') 

disp(Bmat) 

disp{' ') 

disp('press any key to continue...') 

pause 

clc 

elseif choice==4, % eigenvalues and root loci 

clc 

disp(' ') 

disp('After you view the root loci plot, a meta file is made.') 
disp('When you are done a screen will tell you the file names') 
disp('of the meta files. To get a hard copy of the plots, you') 
disp('must graphics post process (GPP) the files for your') 
disp('particular printer set-up then, print.') 
disp(' ') 

disp('NOTE: If ALL roots are real. MATLAB will NOT plot them') 
disp(' in the A->- tand plane, but will plot the root against') 

disp(' its position in the vector (e.g. the first root would') 

disp(' be ploted as (l.root))') 

disp (' ') 

disp('press any key to continue...') 

pause 

clc 

disp(' ') 

disp(' *** EIGENVALUES ***') 

eval(['disp('' ',filenamel,''')']) 

disp(' ') 

disp('Uncoupled') 

disp(’ ') 

disp('Longitudinal plant') 
disp(' ') 

disp(' Root wn damping') 

[wn,z]zdamp(Rlonaug),- 
disp([Rlonaug,wn,z]) 
disp(' ') 

disp('Lateral/Directional plant') 

[wn, z] =damp (Rlataug) ; 
disp(' ') 

disp(' Root wn damping') 

disp([Rlataug,wn,z]) 
disp (' ') 

disp('press any key to continue...') 
pause 

plot(Rlonaug,'*').grid,title{'Roots of Longitudinal Plant') 
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pause 

!del rootlon.* 
meta rootlon 

plot(Rlataug,'*').grid,title('Roots of Lateral/Directional Plant') 

pause 

!del rootlat.* 
meta rootlat 
clc 

disp(' *** EIGENVALUES ***') 

eval(['disp('' filenamel) 

disp(' ') 

disp('Coupled Plant’) 
disp(' ') 

[wn, z] xdamp (Rcoup) ; 

disp(* Root wn damping') 

disp([Rcoup wn z]) 
disp(' ') 

disp('press any key to continue...') 

pause 

plot(Rcoup,'*').grid,title('Roots of Coupled Plant') 
pause 

!del rootcoup.* 
meta rootcoup 
clc 

disp(' ') 

disp ('Plots are saved under the following filenames:') 
disp(' ') 

disp ('Longitudinal roots • rootlon.met') 
disp('Lateral/Directional roots - rootlat.met') 
disp('Coupled roots - rootcoup.met') 
disp(' ') 

disp('press any key to continue...') 
pause 

elseif choice==5, % key control parameters 

clc 

disp(' ') 

disp(' *** KEY CONTROL PARAMETERS (screen 1 of 2) ***') 

eval(['disp('' ',filenamel,''')']) 

disp(' ') 
fprintf(' 
disp(' ') 
disp(' 
fprintf(' 
fprintf(' 
fprintf(' 
disp(' ') 
disp(' 
fprintf(' 
fprintf(' 
fprintf(' 
disp(' ') 
disp(' ') 

disp('press any key to continue...') 
pause 
clc 

disp(' ') 

disp(' *** KEY CONTROL PARAMETERS (screen 2 of 2) ***') 

eval ([' disp (' ' ', filenamel,'")']) 

disp(' ') 

disp(' Cooper Harper Pilot Ratings') 


ss coupling * %6.2f \n',xcouple) 


Designed damping') 

pitch * %6.1f ft-lbs/(rad/sec)\n',desdmdq) 
roll = %6.If ft-lbs/(rad/sec)\n',desdrdp) 
yaw * %6 .If ft-lbs/(rad/sec)\n',desdndr) 

Control Power') 

pitch * %6.1f ft-lbs/in\n',cppitch) 
roll = %6.If ft-lbs/in\n',cproll) 
yaw > %6.If ft-lbs/in\n',cpyaw) 
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disp(' damping/moment of inertia') 

fprintf('pitch (dM/dq)/lyy » %6.2f [ft-lbs/(rad/sec)]/(slug 

ft A 2)\n'.prpitch) 

fprintf(' roll (dR/dp)/Ixx » %6.2f [ft-lbs/(rad/sec)]/(slug 

ft A 2)\n',prroll) 

fprintf(' yaw (dN/dr)/Izz * %6.2f [ft-lbs/(rad/sec)]/(slug 

f t A 2)\n',pryaw) 
disp(' ') 

disp(' control power/moment of inertia') 

fprintf('pitch (dM/in)/Iyy * %6.2f (ft-lbs/in)/(slug ft A 2)\n',cpipitch) 
fprintf(' roll (dR/in)/Ixx - %6.2f (ft-lbs/in)/(slug ft A 2)\n',cpiroll) 
fprintf(' yaw (dN/in)/Izz * %6.2f (ft-lbs/in)/(slug ft A 2)\n',cpiyaw) 
disp(' ') 
disp(' ') 

disp('press any key to continue...') 
pause 

elseif choicest, % command bandwidth plots 

clc 

disp(' ') 

disp('After you view a bode plot of the transfer function from') 
disp('input to state output, a meta file is made. When you exit,') 
disp('a screen will tell you the file names of the meta files.') 
disp('To get a hard copy of the plots, you must graphics post') 
disp('process (GPP) the files for your particular printer set-up') 
disp('then print.') 
disp(' ') 

if Vinf<20 
cmdbwplh 
else, 

cmdbwplc 

end 

elseif choice==0, 
flag=0; 
else 

disp (' ') 

disp('Enter a number on the menu') 
pause(3) 
end 
end 
end 
% 

% *** output to disk (text file) *** 

% 


diary on 
diary off 
delete diary 
diary on 
disp(' 

eval (['disp('' 
disp(' ') 
disp(' 
disp(' ') 
disp(' 
disp(' ') 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
disp(' ') 
disp(' 


*** RESULTS ***') 

',filenamel,''')']) 

*** INPUT DATA ***') 

Flight Conditions') 

Forward velocity = %6.0f kts\n',Vinf/I.69) 
Temperature ■ %6.0f degs F\n',temp) 
Pressure altitude * %6.0f ft\n',PA) 
Auxiliary thrust « %6.0f lbs\n',Taux) 

Fuselage') 
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disp (' ') 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
disp(' ') 
disp(' 
disp(' ') 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
if airfoil 
disp(' 
else 

disp(' 

end 


Gross weight » %6.0f lbs\n',GW) 
Equivalent flat plate area » %6.lf ft A 2\n',Afh) 
Vertical projected area * %6.1f ft A 2\n',Afv) 

CG height above waterline » %6.1f ft\n',zcg) 

CG fuselage station * %6.1f ft\n',xcg) 

CG position rt of buttline - %6.1f ft\n',ycg) 

Ixx ■ %6. if slug £t A 2\n*, Ixx) 
lyy « %6. If slug ft A 2\n',Iyy) 
Izz - %6. If slug f t A 2\n', Izz) 
Ixz * %6. If slug ft A 2\n',Ixz) 
Downwash ratio - %6.2f \n',vfvl) 

Main Rotor') 


-- 1 . 


Number of blades * %6.0f \n',b) 

Rotor radius * %6.1f ft\n',R) 

Blade chord = %6.1f ft\n',cblade) 

Blade twist ■ %6.2f degs\n',twist*57.3) 

Blade airfoil = HH-02') 


Blade airfoil = VR-12') 


fprintf(' Blade lift curve slope = %6.2f \n',a) 

fprintf(' Blade weight ■ %6.1f lbs\n',wblade) 

fprintf(' Rotational velocity « %6.2f rads/sec\n',omega) 

fprintf(' Blade grip length = %6.1f ft\n',grip) 

fprintf(' Hinge offset * %6.1f ft\n',e) 

fprintf(' Flapping moment of inertia ^ %6.1f slug ft A 2\n',Ib) 

fprintf(' Hub height above waterline = %6.1f ft\n',hmd) 

fprintf(' Hub fuselage station * %6.1f ft\n',lmd) 

fprintf(' Hub position rt of buttline « %6.1f ft\n',ymd) 
fprintf(' Mast incidence « %6.2f deg\n',im*57.3) 

disp(' ') 
if notar==0 


disp(' Tail rotor (zero if NOTAR)') 

disp(' ') 

fprintf(' Number of blades * %6.1f \n',bt); 

fprintf(' Blade chord * %6.1f ft\n',cot) 

fprintf(' Blade radius = %6.1f ft\n',Rt) 

fprintf(' Lift curve slope = %6.2f \n',at) 

fprintf(' Rotational velocity * %6.2f rad/sec\n',ohmt) 

fprintf (' Flapping moment of inertia * %6.1f slug f t'‘2\n', Ibt) 
fprintf(' Delta-3 angle = %6.2f deg\n',delta3*57.3) 

fprintf(' Blade twist * %6.2f deg\n',thetalt*57.3); 

fprintf(' Hub height above waterline * %6.1f ft\n',htd) 
fprintf(' Hub fuselage station > %6.1f ft\n',ltd) 

fprintf(' Hub position rt of buttline * %6.1f ft\n',ytd) 
disp(' ') 
elseif notar»l 


disp(' NOTAR') 

disp(' ') 

fprintf(' Height edsove waterline * %6.1f ft A 2\n',htnd) 

fprintf(' Fuselage station ■ %6.1f ft A 2\n',ltnd) 

fprintf(' Position right of buttline « %6.1f ft A 2\n',ytnd) 
fprintf(' NOTAR boom diameter = %6.lf ft A 2\n',dian) 

fprintf(' Swirl angle at boon * %6.2f deg\n',swirl*57.3) 

fprintf(' Maximum thruster force » %6.1f lbs\n',Ytmaxn) 

fprintf(' Thrust fuselage station = %6.1f ft A 2\n',lttnd) 
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dispt' ' 
end 
diep(' 
disp(' ') 
fprintf{' 
fprintf{' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fpririwf (' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
disp(' ') 
disp (' 
disp(' ') 
fprintf (' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
disp(' ') 
disp(' 
disp(' ') 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
fprintf(' 
disp(' ') 
disp(' 
disp(' ') 
fprintf(' 
fprintf(' 
fprintf(' 


Wing') 

Area ■ 
Span = 
CL m 
CDO m 

Tip cord * 
Root cord «= 
Wing efficiency factor » 
Zero lift angle « 
Angle of incidence • 
Lift curve slope « 
Height above waterline ■ 
Fuselage station » 
Position right of buttline * 
Rotor downwash ratio * 
Fuselage downwash ratio « 


%6 . If ft'^Xn' , Swing) 

%6.If ft\n',bwing) 

%6.2f \n',CLwing) 

%6.4f \n’,CDowing) 

%6.1f ft\n',ctw) 

%6 .If ft\n',crw) 

%6.2f \n',ewing) 

%6.2f deg\n',alplow*57.3) 
%6.2f deg\n',iw*57.3) 

%6.2f \n',aw) 

%6.if ft\n',hwd) 

%6.If ft\n',lwd) 

%6.If ft\n*,ywd) 

%6.2f \n',vwvl) 

%6.2f \n',detafdalpfw) 


Horizontal tail' 


Area 

Span 

CL 

CDo 

Zero lift angle 
Angle of incidence 
Lift curve slope 
Height above waterline 
Fuselage station 
Position right of buttline 
Dynamic pressure ratio 
Rotor downwash ratio 
Fuselage downwash ratio 


* %6.1f ft^Xn' ,Shoriz) 

* %6.If ft\n',bhoriz) 

* %6.2f \n',CLhoriz) 

* %6.4f \n',CDohoriz) 

« %6.2f deg\n',alploh*57.3) 

* %6.2f deg\n',iw*57.3) 

= %6.2f \n',ah) 

■ %6.If ft\n',hhd) 

= %6.If ft\n',lhd) 

« V6.If ft\n',yhd) 

* %6.2f \n',qhq) 

= »6.2f \n',vhvl) 

* %6.2f \n',detafdalpfh) 


Vertical tail' 


Area 

Span 

CL 

CDo 

Height above waterline 
Fuselage station 
Position right of buttline 
Zero lift angle 
Maximum Cl 
Dynamic pressure ratio 
Lift curve slope 


%6.1f ff'ZXn',Svert) 

%6.1f ft\n',bvert) 

%6.2f \n',CLvert) 

%6.4f \n',CDovert) 

%6. If ft\n',hvd) 

%6. If ft\n',lvd) 

%6. If ft\n',yvd) 

%6.2f deg\n',alplov*57.3) 
%6.2f \n*,clvertmax) 

%6.2f \n',qvq) 

%S.2f \n',av) 


Rigging') 


fprintf(' Long cyclic pitch/inch defl * %6.2f deg/in\n',dblmddele*57.3) 
fprintf(’ Lat cyclic pitch/inch defl * %6.2f deg/in\n',dalmddela*57.3) 
fprintf (' Collective pitch/inch defl * %6.2f 

deg/in\n',dthetomddelc*57.3) 

fprintf (' Tail rotor pitch change/defl * %6.2f 

deg/unit\n',dthetotddelp*57.3) 
disp(' Max deflection of control') 

fprintf(' from neutral for NOTAR * %6.2f units\n',dphinddelp*57.3) 


pitch change/defl 
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disp(' Displacement of anti - torque') 

fprintf(' control until full rudder * %6.2f units\n•,maxr) 
disp(' ') 

disp(' *** CALCULATED DATA ***') 

disp(' ') 

disp(' Main Rotor') 

disp(' ') 

fprintf(' thrust * %6.1f lbs\n',T) 

fprintf(' torque « %6.1f ft-lbs\n'.Qrotor) 

fprintf(' advance ratio * %6.lf \n',mu) 

fprintf(' inflow parameter wrt TPP ■ %6.3f \n',lamp) 

fprintf(' Tip path angle * %6.1f degs\n',altpp*57.3) 

fprintf(' Rotor coning angle * %6.1f degs\n',ao*57.3) 

fprintf!' 1st lat cyclic tenn-Al « %6.2f degs\n',A1*57.3) 

fprintf!' 1st long cyclic term-Bl * %6.2f degs\n',B1*57.3) 

fprintf!' lateral flapping « %6.1f degs\n',bls*57.3) 

fprintf(' longitudinal flapping * %6.lf degs\n',als*57.3) 

fprintf(' Lock number * %6.lf \n',lockno) 

disp(' ') 

disp(' Tail Rotor {zero if NOTAR)') 

disp(' ') 

fprintf(' tail rotor thrust * %6.1f lbs\n',Tt) 

fprintf(' advance ratio * %6.1f \n',mut) 

fprintf(' inflow parameter * %6.3f \n',lampt) 

fprintf(' Rotor coning angle * %6.1f degs\n’,aot*57.3) 

fprintf(' lateral flapping « %6.2f degs\n',blst*57.3) 

fprintf(' longitudinal flapping = %6.2f degs\n',alst*57.3) 

fprintf(' Lock number = %6.1f \n',locknot) 

disp(' ') 
disp(' ') 

disp(' State Matrices') 

disp(' ') 

disp('Longitudinal uncoupled plant (A or F depending on notation)') 

disp('States are [u w q theta]') 

disp(' ') 

disp(Flonaug) 

disp(' ') 

disp(' ') 

disp ('Longitudinal uncoupled input matrix (B or G depending on notation) ') 

disp('Inputs are [longitudinal cyclic, collective, lateral cyclic, 

pedals]') 

disp(' ') 

disp(Glonaug) 

disp(' ') 

disp(' ') 

disp (' Lateral/directional uncoupled plant (A or F depending on notation) ') 

disp('States are [v p phi r psi]') 

disp(' ') 

disp(Flataug) 

disp(' ') 

disp(' *) 

disp('Lateral/directional uncoupled input matrix (B or G depending') 
disp('on notation)') 

disp('Inputs are [longitudinal cyclic, collective, lateral cyclic, 

pedals]') 

disp(Glataug) 

disp(' ') 

disp(' ') 

disp('Coupled plant (A or F depending on notation)') 
disp('States are [u w q theta v p phi r psi]') 
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disp(' ') 
disp(Amat) 
disp(' ') 
disp(' ') 

disp('Coupled input matrix (B or G depending on notation)') 

disp('Inputs are [longitudinal cyclic, collective, lateral cyclic, 

pedals]') 

disp(' ') 

disp(Bmat) 

disp(' ') 

disp(' ') 

disp(' Eigenvalue') 

disp(' ') 

disp('Uncoupled') 

disp(' ') 

disp('Longitudinal plant') 
disp (' ') 

disp(' Root wn damping') 

[wn, z] »darrtp (Rlonaug) ; 
disp ([Rlonaug,wn,z]) 
disp(' ') 

disp('Lateral/Directional plant') 

[wn,z]xdamp(Rlataug); 
disp(' ') 

disp(‘ Root wn damping') 

disp([Rlataug,wn,z]) 
disp(' ') 
disp (' ') 

disp('Coupled Plant') 
disp(' ') 

[wn,z]sdamp(Rcoup); 

disp(' Root wn damping') 

disp([Rcoup wn z]) 
disp (' ') 
disp (' ') 

disp(' *** KEY CONTROL PARAMETERS ***') 

disp(' ') 

fprintf(' cross coupling = %6.2f \n',xcouple) 
disp(' ') 

disp(' Designed dancing') 

fprintf(' pitch = %6.1f ft-lbs/(rad/sec)\n',desdmdq) 

fprintf(' roll = %6.1f ft-lbs/(rad/sec)\n',desdrdp) 

fprintf(' yaw = %6.1f ft-lbs/(rad/sec)\n'.desdndr) 

disp(' ') 

disp(' Control Power') 

fprintf(' pitch = V6.1f ft-lbs/in\n',cppitch) 

fprintf{' roll = %6.1f ft-lbs/in\n',cproll) 

fprintf(' yaw = %6.1f ft-lbs/in\n',cpyaw) 

disp(' ') 

disp(' Cooper Harper Pilot Ratings') 

disp(' damping/moment of inertia') 

fprintf('pitch (dM/dq)/Iyy = %6.2f [ft-lbs/(rad/sec)]/(slug 

ft A 2)\n',prpitch) 

fprintf(’ roll (dR/dp)/Ixx = %6.2f [ft-lbs/(rad/sec)]/(slug 

ft A 2)\n',prroll) 

fprintf(' yaw (dN/dr)/Izz = %6.2f [ft-lbs/(rad/sec)]/(slug 

ft A 2)\n',pryaw) 
disp(' ') 

disp(' control power/moment of inertia') 

fprintf('pitch (dM/in)/Iyy = %6.2f (ft-lbs/in)/(slug ft A 2)\n',cpipitch) 
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fprintf(' roll (dR/in)/Ixx « %6.2f (ft-lbs/in)/(slug ft“2)\n',cpiroll) 
fprintf(' yaw (dN/in)/Irz * %6.2f (ft-lbs/in)/(slug ft“2)\n',cpiyaw) 
disp(' ') 
disp{' ') 
diary off 
% 

clc 

disp (' ') 

disp(' *** OUTPUT DATA INSTRUCTIONS (screen 1 of 3) ***') 

disp(' ') 

disp('Because this subroutine generates a large number of single') 
disp('value data not shown on the output screen, a text file') 
disp('VARLIST.TXT is on this disk which lists the variable names') 
disp('for all the stability derivatives. Stability derivative') 
disp('contributions for all major aircraft components can be found') 
disp('by reading the text file VARLIST.TXT, then asking MATLAB the') 
disp('variable name corresponding to the derivative.') 
disp(' ') 

disp('Press any key to continue') 
pause 
clc 

dispC ') 

disp(' *** OUTPUT DATA INSTRUCTIONS (screen 2 of 3) ***') 

disp(' ') 

disp('A. Data from the output screen saved to a file named:') 
eval(['disp('' ',filenamel,'.stb'')']) 

disp(' ') 

disp('B. This is a text file, use the TYPE command to view the file') 
disp(' or use a text editor to view/print the file.') 
eval(['flag=exist(''',filenamel,'.stb'');']); 
if flag < 1, 

eval(['!rename diary ',filenamel,'.stb''']); 
else 

eval(['!del ',filenamel,'.stb''']); 
eval(['!rename diary ',filenamelstb'''J); 
end 

disp(' ') 

disp('C. Matrix and vector data saved to a default file named:') 
disp(' mstabdat.mat') 
dispC ') 

disp('D. This is a ".mat" binary file, use the LOAD command to') 
dispC retrieve the data for plotting.') 
disp(' ') 

dispCE. Rename "mstabdat.mat" to another ".mat" file.') 
disp(' The file "mstabdat.mat" will be overwritten when') 
disp(’ the program is executed.') 
disp(' ') 

eval(['disp(''F. Do not rename the file as "',filenamel,'.mat"'')']) 
eval(t'dispt'' The file filenamel,'.mat" is already on disk'')']) 
disp(' and used for future editing.') 
disp(' ') 

disp('Press any key to continue') 
pause 

% 

% *** Output to disk (.mat file containing matrix variables 

% Amat Bmat Rcoup Flataug Glataug Rlataug Plataug Flataug 
% Glataug Rlonaug Plonaug 

% 

% *** Configuring variables for output *** 

% 











save mstabdat Amat Bmat Rcoup Flataug Glataug Rlataug Plataug Flataug ... 
Glataug RXonaug Plonaug 

% 

clc 

disp(' ') 

diap(' *** OUTPUT DATA INSTRUCTIONS (screen 3 of 3) ***') 

disp(' ') 

disp('A. Single value data saved to a default file named:') 
disp(' vstabdat.mat') 
disp(' ') 

disp('B. This is a ".mat" binary file, use the LOAD command to') 
disp(' retrieve the data for plotting.') 
disp(' ') 

dispCC. Rename "vstabdat.mat" to another ".mat" file.') 
disp(' The file "vstabdat.mat" will be overwritten when') 
disp(' the program is executed.') 
disp(' ') 

eval(['disp(''D. Do not rename the file as "',filenamel,'.mat*'')']) 
eval(['disp('' The file "',filenamel,'.mat" is already on disk'')']) 
disp(' and used for future editing.') 

t 

% *** Configuring variables for output *** 

% 

clear Amat Bmat Rcoup Flataug Glataug Rlataug Plataug Flataug ... 
Glataug Plonaug Plonaug num den vAl vBl vals vbls vmu vtheta7 ... 
vao wl vlamp vthetao vctsig vcqsig vchsig valtpp 

% 

save vstabdat 

clear 

clc 

disp(' ') 

disp ( ' *** END STABILITY AND CONTROL ROUTINE ***') 

disp(' ') 

disp('press any key to continue...') 
pause 

format loose 

% 

% return to JANRAD.M 
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N. CMDBWPLH.M 
% CMDBWPLH.M 

% open loop response plots for longitudinal and lateral plants % 
disp('While viewing a plot, press any key to go to the next plot') 
disp <' ') 

disp('Do you want to see longitudinal or lateral/directional plots?') 
disp(' ') 

disp('l. Longitudinal (eight plots total).') 
disp('2. Lateral Directional (ten plots total).') 
disp(' ’) 

pview»input('Enter a number : '); 
clc 

wlogspace(-2,2); 
if pview==l 

% open loop response plots for longitudinal and lateral plants 
w*logspace(-2,2); 

Du=[0 C 0 0] ; 

Cu= [1 Ci 0 0] ; 

Cthet=l0 0 0 1]; 

Cqrat=[0 0 10]; 

Cw=[0 10 0] ; 

disp('longitudinal cylic') 

% command bw e to u 

[NUM,DEN]=ss2tf(Flonaug,Glonaug,Cu,Du,1); 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

title('Open loop response Longitudinal Cyclic to U, Hover') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

!del cbe2uh.* 
meta cbe2uh 

% command bw e to theta 

[NUM,DEN]«ss2tf(Flonaug,Glonaug,Cthet,Du,1); 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

title('Open loop response Longitudinal Cyclic to Theta, Hover') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

!del cbe2theh.* 
meta cbe2theh 
% command bw e to q 

[NUM,DEN] =ss2tf (Flonaug,Glonaug, Cqrat,Du, 1) ,- 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

title('Open loop response Longitudinal Cyclic to Pitch Rate, Hover') 

xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 

pause 

!del cbe2qh.» 
meta cbe2qh 
% command bw e to w 

[NUM,DEN]=ss2tf(Flonaug,Glonaug,Cw,Du,1); 

semilogx(w,20*logl0(bode(NUM,DEN, w))) , grid 

title('Open loop response Longitudinal Cyclic to W, Hover') 

xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 

pause 

(del cbe2wh.* 
meta cbe2wh 

% 

% now collective 
disp('collective') 

% command bw c to u 

[NUM,DEN]=ss2tf(Flonaug,Glonaug,Cu,Du, 2); 
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semilogx(w,20*logl0(bode(NUM, DEN, w))} .grid 
title('Open loop response Collective to U, Hover') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

!del cbc2uh.* 
met a cbc2uh 

% command bw c to theta 

[NUM,DEN]»ss2tf(Flonaug.Glonaug,Cthet.Du,2); 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 
title('Open loop response Collective to Pitch, Hover') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

!del cbc2theh.* 
meta cbc2theh 
% command bw c to q 

[NUM,DEN]=ss2tf(Flonaug,Glonaug,Cqrat,Du,2); 

semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

title('Open loop response Collective to Pitch Rate, Hover') 

xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 

pause 

!del cbc2qh.* 
meta cbc2qh 
% command bw c to w 

[NUM,DEN]=ss2tf(Flonaug.Glonaug,Cw,Du,2); 

semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

title('Open loop response Collective to W, Hover') 

xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 

pause 

Jdel cbc2wh.* 
meta cbc2wh 
clc 

disp(' ') 

disp('Plots are Baved under the following filenames:') 
disp(' ') 

disp('Longitudinal Cyclic') 

disp('Longitudinal Cyclic to U, Hover - cbe2uh.met') 

disp('Longitudinal Cyclic to Theta, Hover - cbe2theh.met') 

disp('Longitudinal Cyclic to Pitch Rate, Hover - cbe2qh.met') 

disp('Longitudinal Cyclic to W, Hover - cbe2wh.met') 

disp ('') 

disp('Collective') 

disp('Collective to U, Hover - cbc2uh.met') 
disp('Collective to Pitch, Hover - cbc2theh.met') 
disp('Collective to Pitch Rate, Hover - cbc2qh.met') 
disp('Collective to W, Hover - cbc2wh.met') 
disp(' ') 

disp('press amy key to continue ...') 
pause 

% 

% now for lateral directional plant 
elseif pview==2 

% 

% now for lateral directional plaint 
Cphi= [0 0 1 0 0]; 

Cv= [1 0 0 0 0]; 

Cp= [0 1 0 0 0]; 

Cr=[0 0 0 1 0];% yaw rate 
Cy= [0 0 0 0 1] ;% yaw auigle 
% lateral cyclic 
disp('lateral cyclic') 





% command bw lateral cyclic to bank 

[NUM,DEN]»ss2tf(Flataug,Glataug,Cphi,Du,3); 

semilogx(w,20*logl0(bode(NUM,DEN,w))) , grid 

title('Open loop response Lateral Cyclic to Bank, Hover') 

xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 

pause 

!del cba2phih.* 
meta cba2phih 

% command bw lateral cyclic to sideslip (lateral velocity) 

[NUM,DEN]=»ss2tf(Flataug,Glataug,Cv,Du, 3); 
semilogx(w,20*logl0(bode(NDM,DEN,wj)).grid 

title('Open loop response Lateral Cyclic to Sideslip (v), Hover') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

(del cba2vh.* 
meta cba2vh 

% command bw lateral cyclic to roll rate 
[NUM,DEN]sss2tf(Flataug,Glataug,Cp,Du,3); 
semilogx(w,20*logl0(bode(NUM,DEN,w)}).grid 

title('Open loop response Lateral Cyclic to Roll Rate, Hover') 
xlabel('Angular Frequency (rad/sec)').ylabel('Gain (dB)') 
pause 

!del cba2ph.* 
meta cba2ph 

% command bw lateral cyclic to yaw rate 
[NUM,DEN]=ss2tf(Flataug,Glataug,Cr,Du,3); 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

title(’Open loop response Lateral Cyclic to Yaw Rate, Hover') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

(del cba2rh.* 
meta cba2rh 

% command bw lateral cyclic to yaw angle 

[NUM,DEN]=ss2tf(Flataug,Glataug,Cy,Du,3); 

semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

title('Open loop response Lateral Cyclic to Yaw, Hover') 

xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 

pause 

!del cba2yh.* 
meta cba2yh 

% 

% pedals 

% command bw pedals to bank 
disp('pedals') 

[NUM,DEN]=ss2tf(Flataug,Glataug(:,4)*30,Cphi, [0],1); 

semilogx(w,20*logl0(bode(NUM,DEN,v))),grid 

title('Open loop response Pedals _o Bank, Hover') 

xlabel('Angular Frequency (rad/sec)').ylabel('Gain (dB)') 

pause 

(del cbp2phih.* 
meta cbp2phih 

% command bw pedals to sideslip 

[NUM.DEN]=ss2tf(Flataug,Glataug(:,4)*30,Cv,[0],1); 

semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

title('Open loop response Pedals to Sideslip (v), Hover') 

xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 

pause 

!del cbp2vh.* 
meta cbp2vh 

% command bw pedals to roll rate 
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[NUM,DEN]x8s2t£(Flataug,Glataug(:,4)*30,Cp, [0],1) ; 
semilogx(w,20*logl0(bode(NUM,DEN,w))) ,grid 
title('Open loop response Pedals to Roll Rate, Hover') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

!del cbp2ph.* 
meta cbp2ph 

% command bw pedals to yaw rate 

[NUM,DEN] =ss2tf (Flataug,Glataug (: ,4) *30, Cr, [0] , 1) ,- 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 
title('Open loop response Pedals to Yaw Rate, Hover') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

!del cbp2rh.* 
meta cbp2rh 

% command bw pedals to yaw 

[NUM,DEN]ess2tf(Flataug,Glataug(:,4)*30,Cy,[0],1); 

semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

title('Open loop response Pedals to Yaw, Hover') 

xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 

pause 

•del cbp2yh.* 
meta cbp2yh 
clc 

disp(' ') 

disp('Plots are saved under the following filenames:') 
disp(' ') 

disp('Lateral cyclic') 

disp('Lateral Cyclic to Bank, Hover - cba2phih.met') 
disp('Lateral Cyclic to Sideslip (v), Hover - cba2vh.met') 
disp('Lateral Cyclic to Roll Rate, Hover - cba2ph.met') 
disp('Lateral Cyclic to Yaw Rate, Hover - cba2rh.met') 
disp('Lateral Cyclic to Yaw, Hover - cba2yh.met') 
disp('Pedals') 

disp('Pedals to Bank, Hover - cbp2phih.met') 
disp('Pedals to Sideslip (v), Hover - cbp2vh.met') 
disp('Pedals to Roll Rate, Hover - cbp2ph.met') 
disp('Pedals to Yaw Rate, Hover - cbp2rh.met') 
disp('Pedals to Yaw, Hover - cbp2yh.met') 
disp(' ') 

disp('press any key to continue ...') 
pause 

% 

else 

end 
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N. 


CMDBWPLC .M 


% CMDBWPLC.M 

% open loop response plots for longitudinal and lateral plants 

% 

disp('While viewing a plot, press any key to go to the next plot') 
disp(' ') 

disp('Do you want to see longitudinal or lateral/directional plots?') 
disp(' ') 

disp('l. Longitudinal (eight plots total).') 

disp('2. Lateral Directional (ten plots total).') 

disp(' ') 

pviewainput ('Enter a number : ') ,• 
clc 

wslogspace(-2,2); 
if pview«l 
Du= [0 0 0 0] ; 

Cu= [1 0 0 0]; 

Cthet= [0 0 0 1] ; 

Cqrat*[0010],- 
Cw*[0 10 0] ; 

disp {'longitudinal cylic') 

% command bw e to u 

[NUM,DEN]=ss2tf(Flonaug,Glonaug,Cu,Du,1); 
semilogx(w,20*logl0(bode(NUM,DEN,w))) ,grid 

title('Open loop response Longitudinal Cyclic to U, Cruise') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

!del cbe2uc.* 
meta cbe2uc 

% command bw e to theta 

[NUM,DEN]=ss2tf(Flonaug,Glonaug,Cthet,Du,1); 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

title('Open loop response Longitudinal Cyclic to Theta, Cruise') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

!del cbe2thec.* 
meta cbe2thec 
% command bw e to q 

[NUM,DEN]=ss2tf(Flonaug,Glonaug,Cqrat,Du,1); 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

title('Open loop response Longitudinal Cyclic to Pitch Rate, Cruise') 

xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 

pause 

!del cbe2qc.* 
meta cbe2qc 
% command bw e to w 

[NUM,DEN]sss2tf{Flonaug,Glonaug,Cw,Du,1); 

semilogx(w,20*logl0(bode(NUM,DEN,w))).grid 

title('Open loop response Longitudinal Cyclic to W, Cruise') 

xlabel('Angular Frequency (rad/Bec)'),ylabel('Gain (dB)') 

pause 

!del cbe2wc.* 
meta cbe2wc 

% 

% now collective 
disp('collective') 

% command bw c to u 

[NUM,DEN]=ss2tf(Flonaug,Glonaug,Cu,Du,2); 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 








title ('Open loop response Collective to U, Cruiee') 
xlabel('Angular Frequency (red/sec)').ylabel('Gain (dB)') 

pause 

!del cbc2uc.* 
meta cbc2uc 

% command km c to theta 

[NUM,DEN]«ss2tf(Flonaug,Glonaug,Cthet,Du,2); 
semilogx(w,20*logl0(bode(NUM,DEN,w))) ,grid 
title('Open loop response Collective to Pitch, Cruise') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

ldel cbc2thec.* 
meta cbc2thec 
% command bw c to q 

[NUM,DEN]»ss2tf(Flonaug,Glonaug,Cqrat,Du,2); 
semilogx(w,20*logl0(bode(HUM,DEN,w))),grid 

title('Open loop response Collective to Pitch Rate, Cruise') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

!del cbc2qc.* 
meta cbc2qc 
% command bw c to w 

(NUM,DEN]=ss2tf(Flonaug,Glonaug,Cw,Du,2); 

semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

title('Open loop response Collective to W, Cruise') 

xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 

pause 

!del cbc2wc.* 
meta cbc2wc 
clc 

disp(' ') 

disp('Plots are saved under the following filenames:') 
disp(' ') 

disp('Longitudinal Cyclic') 

disp('Longitudinal Cyclic to U, Cruise - cbe2uc.met') 
disp('Longitudinal Cyclic to Theta, Cruise - cbe2thec.met') 
disp('Longitudinal Cyclic to Pitch Rate, Cruise - cbe2qc.met') 
disp('Longitudinal Cyclic to W, Cruise - cbe2wc.met') 
disp('') 

disp('Collective') 

disp('Collective to U, Cruise - cbc2uc.met') 
disp('Collective to Pitch, Cruise - cbc2thee.met') 
disp('Collective to Pitch Rate, Cruise - cbc2qc.met') 
disp('Collective to w, Cruise - cbc2wc.met') 
disp(' ') 

disp('press any key to continue ...') 
pause 
% 

% now for lateral directional plant 
elseif pview«2 
Cphis [0 0 1 0 0]; 

Cv= [1 0 0 0 0]; 

Cp= [0 1 0 0 0]; 

Cr*[0 0 0 1 0];% yaw rate 
Cy=[0 000 1];% yaw angle 
% lateral cyclic 
disp('lateral cyclic') 

% command bw lateral cyclic to bank 

[NUM,DEN]*ss2tf(Flataug,Glataug,Cphi,Du,3); 

semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 




title('Open loop response Lateral Cyclic to Bank, Cruise') 
xlabel('Angular Frequency (rad/sec)').ylabel('Gain (dB)') 

pause 

idel cba2phic.* 
meta cba2phic 

% command bw lateral cyclic to sideslip (lateral velocity) 

[MUM,DBM]-ss2tf(Flataug,Glataug,Cv,Du,3); 

semilogx(w,20*logl0(bode(NUM,DEN,w))) ,grid 

xlabel('Angular Frequency (rad/sec)').ylabel('Gain (dB)') 

titleCOpen loop response Lateral Cyclic to Sideslip (v) , Cruise') 

pause 

idel cba2vc.* 
meta cba2vc 

% command bw lateral cyclic to roll rate 
[MUM,DEN]B8s2t£(Flataug,Glataug,Cp,Du,3); 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

title('Open loop response Lateral Cyclic to Roll Rate, Cruise') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

idel cba2pc.* 
meta cba2pc 

% command bw lateral cyclic to yaw rate 

[MUM,DEN]*ss2tf(Flataug,Glataug,Cr,Du,3); 

semilogx(w,20*logl0(bode(MUM,DEN,w))),grid 

xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 

title('Open loop response Lateral Cyclic to Yaw Rate, Cruise') 

pause 

idel cba2rc.* 
meta cba2rc 

% command bw lateral cyclic to yaw angle 

[NUM,DEN]sss2tf(Flataug,Glataug,Cy,Du,3); 

semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 

xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 

title('Open loop response Lateral Cyclic to Yaw, Cruise') 

pause 

idel cba2yc.* 
meta cba2yc 

% 

% pedals 

% command bw pedals to bank 
disp('pedals') 

[NUM,DEN]sss2tf(Flataug,Glataug(:,4)*30,Cphi, [0],1); 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 
title('Open loop response Pedals to Bank, Cruise') 
xlabel('Angular Frequency (rad/sec)'),ylabel{'Gain (dB)') 
pause 

idel cbp2phic.* 
meta cbp2phic 

% command bw pedals to sideslip 

[NUM,DEN]-ss2tf(Flataug,Glataug(:,4)*30,Cv,[0],1); 

semilogx(w,20*logl0(bode(NUM,DEN,w))) , grid 

title('Open loop response Pedals to Sideslip (v), Cruise') 

xlabel('Angular Frequency (rad/sec)').ylabel('Gain (dB)') 

pause 

idel cbp2vc.* 
meta cbp2vc 

% command bw pedals to roll rate 

[NUM,DEN]*ss2tf(Flataug,Glataug(:,4)*30,Cp,[0],1); 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 
title('Open loop response Pedals to Roll Rate, Cruise') 
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xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 

pause 

!del cbp2pc.* 
meta cbp2pc 

% command bw pedals to yaw rate 

[NUM,DEN]«ss2tf(Flataug,Glataug(:,4)*30,Cr, [0],1); 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 
title('Open loop response Pedals to Yaw Rate, Cruise') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

!del cbp2rc.* 
meta cbp2rc 

% command bw pedals to yaw 

[NUM,DEN]>ss2tf(Flataug,Glataug(:,4)*30,Cy,[0],1); 
semilogx(w,20*logl0(bode(NUM,DEN,w))),grid 
title('Open loop response Pedals to Yaw, Cruise') 
xlabel('Angular Frequency (rad/sec)'),ylabel('Gain (dB)') 
pause 

!del cbp2yc.* 
meta cbp2yc 

% 

clc 

disp(' ') 

disp('Plots are saved under the following filenames:') 
disp(' ') 

disp('Lateral cyclic') 

disp('Lateral Cyclic to Bank, Cruise - cba2phic.met') 
disp('Lateral Cyclic to Sideslip (v), Cruise - cba2vc,met') 
disp('Lateral Cyclic to Roll Rate, Cruise - cba2pc.met') 
disp('Lateral Cyclic to Yaw Rate, Cruise - cba2rc.met') 
disp('Lateral Cyclic to Yaw, Cruise - cba2yc.met') 
disp('Pedals') 

disp('Pedals to Bank, Cruise - cbp2phic.met') 
disp('Pedals to Sideslip (v), Cruise - cbp2vc.met') 
disp('Pedals to Roll Rate, Cruise - cbp2pc.met') 
disp('Pedals to Yaw Rate, Cruise - cbp2rc.met') 
disp('Pedals to Yaw, Cruise - cbp2yc.met') 
disp(' ') 

disp('press any key to continue ...') 
pause 
* 

else 

end 
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0. DKRXV2.M 
% DERIV2.M 

t computes the derivative of a polynomial 

function[deriv]-deriv(polyno); 

der» [2 1 0] ; 

dsder.*polyno; 

deriv-fd(l,1) d(l,2)]; 
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